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FINAL  REPORT  -  ONR  CONTRACT  #N00014- 84-K-0027 

CREEP  AND  FRACTURE  CHARACTERISTICS  OF  MATERIALS  AND  STRUCTURES 
AT  ELEVATED  TEMPERATURES 

PRINCIPAL  INVESTIGATOR  -  DEAN  HAROLD  LIEBOWITZ 
TECHNICAL  DIRECTOR  -  ASSOC.  PROF.  E.  THOMAS  MOYER  JR. 

Significant  research  was  performed  under  ONR  Contract  #N00014-84-K-0027 
during  the  period  of  the  contract.  As  listed  in  Appendix  A,  this  work 
resulted  in  eight  refereed  publications  and  four  invited  presentations  at 
International  Conferences.  The  work  also  resulted  in  four  student  theses 
listed  in  Appendix  B.  In  addition,  experimental  progress  was  made  in 
creep  fracture  testing.  As  outlined  below,  this  work  was  has  not  been 
completed  due  to  the  lack  of  continuation  of  the  contract. 

The  work  under  this  contract  was  concentrated  in  three  major  areas:  the 
effect  of  mixed  mode  loading  on  fracture  characteristics,  the  nature  of 
crack  tip  stress,  strain  and  energy  fields  in  ductile  materials  and  the 
nature  of  crack  tip  stress  strain  and  energy  fields  in  materials 
undergoing  rate  dependent  viscoplastic  deformation.  In  each  of  these 
areas,  new  insight  was  obtained  and  better  understanding  of  the 
fundamental  physical  processes  gained.  (  . 

Early  work  on  this  contract  focused  on  mixed  mode  fracture 
characteristics.  Experimental  studies  and  finite  element  modeling 
determined  specimen  characteristics  and  design  modifications  for  a  mode 
two  fracture  specimen.  This  specimen  has  the  unique  capability  of 
testing  from  pure  mode  one  to  pure  mode  two  without  significant  crack 
face  rotation.  This  work  is  documented  in  two  student  thesis  (students 
were  successful  candidates  for  the  Diplome  degree  through  a  joint, 
cooperative  program  between  the  University  of  Stuttgart  and  the  George 
Washington  University).  In  addition,  further  work  developed  a 
computational  procedure  based  on  the  nodal  force  approach  for  the 
determination  of  stress  intensity  factor  distributions  along  arbitrary 
crack  fronts  in  three  dimensions.  This  work  was  presented  at  an  ASTM 
conference  and  will  appear  in  ASTM  STP  #969.  This  work  is  also 
documented  in  a  student  masters  thesis  (all  student  theses  are  listed  in 
Appendix  B) . 

Work  on  ductile  fracture  was  carried  out  for  three  dimensional,  mode  one 
crack  geometries.  The  effect  of  specimen  thickness  and  material 
hardening  characteristics  was  studied.  In  addition  to  useful 
understanding,  the  thickness  range  where  plane  strain  and  plane  stress 
are  valid  assumptions  were  discovered.  Depending  on  the  ductility  of  the 
material,  the  plane  strain  thickness  did  not  correspond  to  the  ASTM 
requirement  due  to  the  assumptions  of  elasticity  employed.  A  modified 
approach  for  determination  of  plane  strain  thickness  was  proposed.  In 
addition  to  the  three  dimensional  studies  involving  a  stationary  mode  one 
crack,  further  research  was  performed  in  the  areas  of  mixed  mode  ductile 
fracture  and  ductile  crack  growth.  The  mixed  mode  ductile  fracture 
studies  demonstrated  the  crack  opening  characteristics  as  a  function  of 
®o',»  ratio.  It  w a_  demonstrated  that  for  dominant  mode  one,  a  distinct 
notch  effect  is  observed.  This  notch  opening  removes  the  HRR  singularity 


and  produces  a  ductile  zone  which  is  characterized  by  a  weaker  energy 
singularity  than  was  previously  known.  In  addition,  as  mode  two  becomes 
more  dominant,  the  deformed  crack  remains  sharper  and  the  local  HRR 
characteristics  return.  In  addition,  the  significant  rotation  occurs  at 
the  crack  tip  altering  the  amplitude  of  the  local  field.  Mode  one  crack 
growth  studies  were  performed.  This  work  demonstrated  a  computational 
approach  for  accurately  modeling  stable  crack  growth  with  a  commercial 
finite  element  code.  The  physical  results  demonstrate  the  disappearance 
of  the  HRR  zone  due  to  notch  opening  and  the  appearance  of  a  significant 
transition  zone  which  dominates  the  local  fracture  zone.  This  zone  is 
characterized  by  an  energy  singularity  which  is  weaker  than  1/r.  This  is 
a  new  result  which  is  under  further  investigation.  The  two  dimensional 
studies  are  documented  in  student  theses  listed  in  Appendix  B.  To  date, 
work  is  continuing  on  this  problem  and  the  results  are  not  yet  available 
in  the  open  literature.  The  student  theses,  therefore,  are  inclv’.ded  as 
Appendices  D  and  E. 

Studies  on  creep  fracture  characteristics  were  the  focus  of  significant 
study  under  this  contract.  Experimental  work  focused  on  crack  growth 
studies  on  IN  718  at  650  degrees  C.  At  this  temperature,  significant 
constituitive  data  was  available.  These  results  demonstrated  that  the  C* 
integral  was  not  employable  as  a  crack  driving  force  measure.  In 
addition,  it  was  determined  that  experimental  scatter  was  due  to  crack 
front  curvature  effects  which  could  be  minimized  through  careful 
experimental  technique.  The  final  results  demonstrated  a  two  stage 
growth  regime  which  was  numerically  fit  to  explicit  time  functions  for 
crack  growth  simulation.  This  work  was  part  of  an  ONR  progress  report 
and  is  included  as  Appendix  E.  Finite  element  studies  of  creep  crack 
growth  were  performed  and  the  results  are  part  of  a  recent  publication 
included  as  Appendix  F.  This  study  demonstrated  the  influence  of  finite 
strains  in  the  crack  region  and  the  inability  of  local  asymptotic 
solutions  to  characterize  the  stress  fields  near  stationary  and  growing 
cracks.  In  addition,  convergence  and  accuracy  of  the  numerical  approach 
was  studied  extensively.  New  understanding  of  convergence 
characteristics  was  obtained. 

Experiments  were  initiated  at  550  degrees  C  to  determine  if  the  results 
at  650  degrees  were  characteristic  of  creep  crack  growth  in  general  or 
were  a  qualitative  function  of  temperature.  Unfortunately,  insufficient 
constituitive  data  was  available  at  550  degrees  and  the  data  was  not  able 
to  be  analyzed.  Constituitive  tests  were  initiated,  however,  the 
contract  resources  were  not  sufficient  to  complete  the  work.  If  future 
funding  is  available  for  this  work,  the  tests  will  be  completed  and  the 
results  will  be  forwarded  to  ONR. 

The  appendices  of  this  work  document  the  significant  research 
contribution  that  was  made  under  ONR  contract  #N00014-84-K-0027 .-  All 
publications  cited  in  the  appendices  were  forwarded  to  ONR  at  publication 
and  were  incLuded  in  the  quarterly  progress  reports.  In  addition,  all 
publications  have  been  sent  to  the  ONR  distribution  list  at  the  time  of 
publication . 
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Materials  Science.  ASM  International,  1988. 

8]  "Effect  of  Specimen  Thickness  on  Crack  Front  Plasticity 
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Society  of  India,  Vol.  36,  pp .  17,  1984. 
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FRACTURE  TESTS  ON  IN  718 


INSIGHT  INTO  CREEP  FRACTURE  PHENOMENA 

E.  T.  Moyer,  Jr.  and  H.  Liebowitz 
The  George  Washington  University 

ABSTRACT 

Fracture  tests  on  IN  718  superalloy  demonstrate  that  the 

* 

C  fracture  parameter  is  not  a  sufficient  quantity  for  the 
quantitative  description  of  creep  crack:  growth.  The  results 

contained  in  this  communication  show  that  the  crack,  velocity  is 

* 

not  uniquely  predicted  by  C  but  is  also  a  function  of  test 
load.  In  addition,  the  results  indicate  that  the  crack, 
velocity  would  also  be  affected  by  geometry  changes  (e.g., 
specimen  size). 

The  results  presented  in  this  communication  also  demon¬ 
strate  that  crack,  growth  initiates  extremely  early  in  the  test 
history.  No  unique  initiation  time  was  identifiable.  Also 
evident  is  a  two  stage  growth  process  with  stage  1  (charac¬ 
terized  by  constant  crack,  velocity)  contributing  significantly 
to  the  total  useful  life  even  at  relatively  high  initial  crack 
velocities  (on  the  order  of  0.001  inches/minute). 

Investigation  was  made  into  the  widely  observed  scatter  in 
creep  fracture  data  reported  in  the  literature.  This  scatter 
is  often  suggested  to  be  due  to  crack  tunneling,  material 
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variability,  etc.  The  results  presented  in  this  work,  show  that 
initial  crack,  front  curvature,  irregular  geometries,  forming 
inconsistencies  (e.g.,  rolling  irregularities)  cause  extreme 
scatter  in  experimental  results.  These  irregularities, 
however,  are  observable  continuum  phenomena  which  are  incon¬ 
sistent  with  the  assumptions  inherent  in  the  analysis  of  the 
test  data.  When  specimens  exhibiting  these  irregularities  are 
removed  from  the  data  base,  scatter  is  reduced  to  acceptable 
levels  (e.g.,  less  than  10%  in  measured  quantities).  The 
fracture  surfaces  also  indicate  that  tunneling  does  not  occur 
for  the  geometry,  loading,  temperature  and  material  conditions 
studied . 
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FRACTURE  TESTS  ON  IN  718 


INSIGHT  INTO  CREEP  FRACTURE  PHENOMENA 

E.  T.  Moyer,  Jr.  and  H.  Liebowitz 
The  George  Washington  University 

A  series  of  constant  load  creep  fracture  tests  were 
performed  on  IN  718  specimens  at  650°C.  The  specimen  geometry 
was  standard  compact  tension  with  dimensions  (a  =  1.0  in., 

W  =  2.0  in.,  B  =  0.4  in.).  Tests  were  run  for  load  levels 
between  1000-1500  lb.  Mouth  opening  displacement  and  crack, 
length  were  monitored  continuously  during  the  test.  To 
establish  a  sharp  initial  crack.,  the  specimens  were  fatigue 

precracked  at  room  temperature  at  15  Ksi  -  /Tn.  Crack  length 
is  measured  optically  to  a  precision  of  0.0015  in. 

The  first  Figure  is  a  plot  of  the  crack  length  vs  time  at 
five  different  load  levels.  All  the  data  indicates  a  two  stage 
crack  growth  process.  The  first  stage  is  characterized  by 
essentially  constant  crack  velocity  (the  linear  portion  of  the 
crack  length  vs  time  curve)  and  the  second  stage  is  character¬ 
ized  by  continuous  acceleration.  For  the  growth  range  studied, 
stage  1  crack  growth  accounts  for  a  significant  portion  of  the 
growth  history  (at  1000  lb.,  stage  1  accounts  for  approximately 
65%  of  the  time  required  to  increase  the  crack  length  40%;  at 
1500  lb.,  stage  1  accounts  for  approximately  40%  of  this  time 
history).  At  the  load  levels  tested,  crack  growth  was  observed 


very  early  in  the  history.  A  unique  "initiation"  time  was  not 
identifiable. 

The  crack,  length  vs  time  curves  clearly  indicate  that 
stage  1  growth  is  evident  even  when  initial  crack,  velocities 
are  of  a  "average"  magnitude.  Previous  studies  have  indicated 
that  stage  1  is  present  only  for  very  slowly  growing  cracks 
[1].  In  the  data  presented  here,  initial  crack,  velocities 
varied  by  almost  an  order  of  magnitude  and  all  tests  exhibit 
stage  1  behavior  for  significant  portions  of  the  growth  history. 

The  second  Figure  is  a  plot  of  crack  velocity  (da/dt)  vs 
C  .  The  formulation  of  Kumar  and  Shih  (K-S)  is  used  for  the 
calculation  of  C  [2].  This  formulation  is  to  be  preferred  to 
the  Harper  and  Ellison  (H-E)  formulation  for  two  reasons: 
first,  the  assumptions  made  in  the  derivation  are  less  restric¬ 
tive  (e.g.,  zone  size  requirements  in  the  H-E  formulation, 
proportioning  of  deformation  due  to  crack  growth  and  creep, 
etc.)  and  second,  because  the  K-S  formulation  requires  only  a 
knowledge  of  the  geometry  and  loading  and  not  the  load  line 
displacement  rate  (which  is  measured  less  accurately).  Indeed, 

for  reasonable  crack  velocities,  the  H-E  formulation  can  be 

* 

shown  to  be  a  measure  of  da/dt  and  not  C  [3]. 

*  * 

The  da/dt  vs  C  plot  demonstrates  that  the  C  is  not  a 

sufficient  parameter  to  describe  crack  growth.  It  has  been 

postulated  in  the  literature  that  da/dt  can  be  uniquely  related 

* 

to  C  ,  independent  of  loading  and  geometry  (see,  for  example, 
[4,5]).  The  results  presented  by  the  authors  demonstrate  that 


3 


v  da/dt  is  a  function  of  the  test  loading  in  addition  to  the 

* 

parameters  involved  in  calculating  C  .  This  is  understand- 

* 

able  as  during  stage  1  growth,  C  is  steadily  increasing  while 

da/dt  remains  constant.  In  addition,  convergence  toward  a 

* 

unique  da/dt  vs  C  relationship  is  not  evident  until  the  c rack, 
velocity  and  length  have  grown  appreciably.  At  the  larger 
da/dt  values,  many  of  the  assumptions  required  for  the 
application  of  the  K-S  formula  become  dubious.  In  the  region 

in  which  the  K-S  assumptions  are  valid,  the  results  demonstrate 

* 

that  C  is  not  a  sufficient  correlating  parameter. 

For  the  range  of  geometry  and  loading  presented  in  this 

work.,  K  is  not  a  viable  fracture  parameter.  At  1000  lb. 

(  loading,  the  stress  intensity  factor  calculated  from  the  mouth 

v 

opening  displacement  reached  a  value  of  27  Ksi  -  /Tn.  prior  to 

crack  growth  where  the  linear  elastic  K  value  was  17  Ksi  -  /Tn. 
corresponding  to  the  load.  It  was  evident,  therefore,  that  the 
creep  deformation  exceeded  the  K  controlled  region  from  the 
start  of  the  test.  K,  therefore  is  not  a  viable  fracture 
parameter  for  this  data. 

Creep  fracture  studies  often  exhibit  much  experimental 
scatter.  Many  reasons  are  proposed  including  environmental 
effects,  deformation  transitions,  tunneling,  mechanism  transi¬ 
tions,  etc.  Data  which  exhibits  large  scatter  cannot  be  used 
%  to  establish  or  reject  the  validity  of  any  theoretical  model  as 

the  error  in  the  data  can  be  on  the  order  of  the  phenomena 
being  described.  To  minimize  scatter  in  our  results,  data  was 
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only  taken  in  the  range  where  oxidation  effects  are  small. 
Oxidation  influences  are  a  function  of  the  test  duration  and 
the  local  strain  state  near  the  crack  [6],  The  presence  of 
oxidation  can  be  seen  on  the  fracture  surface  as  a  change  in 
color  from  the  standard  metallic  color  to  a  blue  color.  Only 
data  obtained  before  the  color  change  became  appreciable  was 
used  in  the  analysis. 

In  addition  to  avoiding  oxidation  driven  growth  data, 
closer  examination  of  the  fracture  specimens  revealed  observ¬ 
able  causes  for  the  scatter  observed.  Several  specimens  which 
exhibited  data  far  from  the  mean  had  extensively  curved  crack 
fronts  after  fatigue  precracking.  These  specimens  tended  to 
exhibit  much  slower  crack  growth  than  those  with  relatively 
straight  crack  fronts.  If  extreme  curvature  was  exhibited  in 
the  fatigue  crack  (greater  than  approximately  1/16  in.)  the 
data  was  rejected. 

Photo  1  and  Photo  2  show  fracture  specimens  whose  data 
were  excepted  (number  1  was  loaded  at  1000  lb.  and  number  2  was 
loaded  at  1500  lb.).  Both  exhibit  typical  fatigue  cracks  which 
produced  consistent  data.  Photo  3  shows  a  specimen  which  was 
loaded  at  1000  lb.  The  data  from  that  test  exhibit  twice  the 
lifetime  of  the  mean  at  that  load.  The  geometric  discontinuity 
introduced  by  a  machining  error  (the  kink  in  the  notch)  caused 
the  crack  to  grow  in  a  uneven  manner  prolonging  life. 

Other  scatter  occurred  due  to  discontinuities  in  the 
material.  Photo  4  shows  a  specimen  which  was  loaded  at 
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1500  lb.  The  fracture  surface  discontinuity  again  prolonged 


the  life  of  the  specimen.  This  discontinuity  is  believed  to  be 


caused  by  imprecision  in  the  rolling  process  during  forming. 


None  of  the  sources  of  scatter  described  could  have  been  pre¬ 


dicted  without  examination  of  the  fracture  surface.  These 


phenomena,  however,  are  not  due  to  material  variability  or 


microstructure.  All  the  observations  are  continuum  irregu¬ 


larities  which  are  inconsistent  with  the  analytical  assumptions 


of  continuum  crack,  growth. 


After  discarding  the  specimens  with  continuum  disconti¬ 


nuities,  the  data  exhibited  very  little  scatter  from  test  to 


test.  The  data  presented  in  this  work,  is  the  average  of  that 


obtained  from  multiple  tests.  The  test  to  test  differences  in 


crack,  length  was  less  than  5%  and  the  difference  in  crack. 


velocities  was  less  than  10%.  Mouth  opening  displacements  were 


contained  within  a  scatter  band  of  approximately  3%  with 


deviations  in  opening  rate  of  approximately  7%.  It  is  felt 


that  these  numbers  accurately  represent  the  "scatter"  which  is 


due  to  testing  configuration,  material  variability  (which 


should  be  small  since  all  specimens  are  from  a  single  batch  of 


material,  were  heat  treated  identically  and  were  cut  in  the 


same  direction  relative  to  the  rolling)  and  microstructure. 


The  data  presented  demonstrate  the  inability  of  either  C 


or  K  to  be  a  valid  constitutive  parameter  for  creep  crack. 


growth  in  IN  718  at  650°C.  In  addition,  these  results  viewed 


with  o^her  investigators'  work  (for  various  materials,  e.g.. 


VlvA  a  AA  A, 


6 


r 


[1,3,6]),  demonstrate  that  a  valid  fracture  parameter  charac¬ 


terizing  creep  crack,  growth  behavior  for  a  realistic  range  of 


geometry  and  loading  has  yet  to  be  found.  In  addition,  new 


insight  into  the  "Sources  of  Scatter"  have  been  identified. 


This  testing  sequence  suggests  that  continuum  reasons  for 


observed  scatter  can  often  be  identified  which  violate  the 


continuity  assumptions  inherent  in  the  test  procedure  and 


analysis . 


•*,  r* 
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In  the  present  work,  local  crack-tip  field  quantities  under  ductile  material  behavior 
were  studied  for  mixed-mode  loading  ranging  from  pure  mode  I  to  pure  mode  II  under  the 
assumption  of  plane  strain.  In  order  to  become  independent  of  a  specific  specimen,  the  local 
crack-tip  region  was  modeled  as  a  disk  with  the  crack  tip  at  its  center.  Based  on  the 
assumption  of  small  scale  yielding,  displacements  evaluated  from  the  linear  elastic  solution 
were  applied  on  the  model  boundary.  For  ten  comparable  cases  of  mixed-mode  loading  the 
body  response  was  calculated  using  the  J2  flow  theory  of  incremental  plasticity  employing 

small  strain  theory.  The  finite  element  mesh  employed  consisted  of  1178  eight  node  plane 
strain  elements  and  3643  nodes. 

In  the  evalution  of  the  results  emphasis  was  placed  on  : 

i  )  The  investigation  of  field  quantities  in  terms  of  their  exposed  singular  behavior, 
magnitude  and  distribution  inside  the  plastic  zone 

i  i )  The  examination  of  the  influence  of  mixed-mode  loading  on  the  singular  behavior  of  the 
field  quantities  and  the  validity  of  the  HRR  singular  field  for  mixed  modes 

i  i  i )  The  discussion  of  the  strain  energy  density  as  a  criterion  predicting  onset  and  direction 
of  crack  growth  for  mixed  mode  loading  with  ductile  material  behavior  and 

iv)  The  determination  of  the  stress  functions  from  the  finite  element  results  and  their 
comparison  with  the  numerical  calculation  of  an  asymptotic  solution. 
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Figure  1  : 
Figure  2  : 


Definition  of  the  three  modes  of  fracture. 


Definition  of  the  crack-tip  stresses,  showing  rectangular  and  polar 
coordinate  components. 

Figure  3  :  Definition  of  crack  angle  and  fracture  angle  in  the  center  cracked  panel 

with  slanted  crack  under  uniaxial  tensile  stress. 

Figure  4  :  K|  -  K||  curve  according  to  the  fracture  criterion  due  to  Sih  [29]. 

Figure  5  :  Idealized  constitutive  material  behavior  of  (i)  incrementally  elastic- 

plastic  material  conforming  to  incremental  theory  of  plasticity,  (ii) 
incrementally  elastic-plastic  material  conforming  to  deformation 
theory  of  plasticity. 

Figure  6  :  Mapping  of  the  eight-node  prabolic  element  from  spatial  coordinates 

(x,  y)  to  local  coordiantes  (£,  n).  Definition  of  the  node  and  integration 
point  numbering. 

Figure  7  :  Contour  Plots  of  the  von  Mises  yield  stress  for  selected  cases  of  mixed¬ 

mode  loading. 

Figure  8  :  Generation  of  the  eight-node  collapsed  element:  (i)  rectangular  eight  - 

node  element,  (ii)  degenerated  triangular  eight-node  element. 

Figure  9  :  Configuration  of  the  degenerated  crack-tip  elements. 

Figure  10  :  Main  fan  of  the  finite  element  model  of  the  specimen,  inner  radius  :  1mm, 
outer  radius  :  10  mm. 

Figure  1 1  :  Boundary  of  the  finite  element  mesh  of  the  specimen,  inner  radius  :  5mm, 
outer  radius  :  100  mm. 

Figure  12  :  Stress  -  strain  curve  of  the  steel  A304. 

Figure  13  :  A  contour  r  around  the  crack  tip  and  parameters  defining  the  J-integral. 

Figure  14  :  Jr  and  ^-integral  values  for  all  mixed-mode  cases  considered. 

Figure  15  :  Illustration  of  the  crack  tip  region  under  small  scale  yielding  condition. 

Figure  16  :  Effective  von  Mises  stresses  along  line  9  =  0°  for  all  mixed-modes  cases 

considered. 


Figure  17  :  Effective  von  Mises  stresses  around  crack  tip.  Radius  r  =  0.4  mm  for 
selected  mixed-mode  cases. 

Figure  18  :  Full  logarithmic  representation  of  effective  von  Mises  stresses  along  line 
0=0°  for  all  mixed-mode  cases  considered. 

Figure  19  :  -  stresses  along  line  0*0°  for  selected  mixed-mode  cases. 

Figure  20  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r-  1  mm. 

Stress  intensity  factor  ratio  :  K|/Kj|  =  2230/0. 

Figure  21  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/K||  =  1927/987. 

Figure  22  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/K(|  =  1683/1252. 

Figure  23  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  Kt/K||  =*  1405/1462. 

Figure  24  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r=»  1  mm. 

Stress  intensity  factor  ratio  :  K|/Kt|  =  772/1774. 

Figure  25  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/Kj|  =  396/1903. 

Figure  26  :  Equivalent  von  Mises  stress  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/K|j  =  0/2018 

Figure  27  :  Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K|/K|(  =  2230/0. 

Figure  28  :  Equivalent  von  Mises  stress  around  crack  tip.  inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K|/K||  =  1927/987. 

Figure  29  :  Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K(/K||  =  1683/1252. 

Figure  30  :  Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K|/K(|  =  1098/1633. 

Figure  31  :  Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  Kj/K||  =  772/1774. 


Figure  32  :  Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm.  Stress  intensity  factor  ratio  :  K|/K||  =  396/1903. 

Figure  33  :  Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm.  Stress  intensity  factor  ratio  :  K|/K||  =  0/2018. 

Figure  34  :  Equivalent  plastic  strain  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/K||  =  2230/0. 

Figure  35  :  Equivalent  plastic  strain  around  crack  tip.  Radius  r»  1  mm. 

Stress  intensity  factor  ratio  :  K|/Kt|  =  1927/987. 

Figure  36  :  Equivalent  plastic  strain  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/K n  -  1098/1633. 

Figure  37  :  Equivalent  plastic  strain  around  crack  tip.  Radius  r=  1  mm. 

Stress  intensity  factor  ratio  :  K|/K||  =  0/2018. 

Figure  38  :  Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K|/KM  ==  2230/0. 

Figure  39  :  Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 

radius  10  mm.  Stress  intensity  factor  ratio  :  K|/KM  =  1927/987. 

Figure  40  :  Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K(/K u  =  1405/1462. 

Figure  41  :  Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  Kj/K u  =  772/1774. 

Figure  42  :  Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm.  Stress  intensity  factor  ratio  :  K|/K n  =  396/1903. 

Figure  43  :  Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 

radius  20  mm.  Stress  intensity  factor  ratio  :  K|/K||  =  0/2018. 

Figure  44  :  Deformed  Mesh,  outer  radius  r=  1mm. 

Stress  intensity  factor  ratio  :  K|/K||  =  2230/0. 

Figure  45  :  Deformed  Mesh,  outer  radius  r=  1mm. 

Stress  intensity  factor  ratio  :  K(/K n  =  1927/987. 

Figure  46  :  Deformed  Mesh,  outer  radius  r=  1mm. 

Stress  intensity  factor  ratio  :  K(/K||  =  1405/1462. 


Figure  47  :  Deformed  Mesh,  outer  radius  r=  1mm. 

Stress  intensity  factor  ratio  :  K(/K||  =  1098/1633. 

Figure  48  :  Deformed  Mesh,  outer  radius  r»  1mm. 

Stress  intensity  factor  ratio  :  Kt/K||  =  0/2018. 

Figure  49  :  Crack  deformation.  Plotted  crack  length  10  mm. 

Stress  intensity  factor  ratio  :  K|/K(|  =  2230/0. 

Figure  50  :  Crack  deformation.  Plotted  crack  length  10  mm. 

Stress  intensity  factor  ratio  :  K|/K||  =  1927/987. 

Figure  51  :  Crack  deformation.  Plotted  crack  length  10  mm. 

Stress  intensity  factor  ratio  :  K|/Kj|  =  1098/1633. 

Figure  52  :  Crack  deformation.  Plotted  crack  length  10  mm. 

Stress  intensity  factor  ratio  :  K|/K||  =  0/2018. 

Figure  53  :  Deformed  and  undeformed  meshes.  Inner  radius  1  mm,  outer  radius  10  mm. 
Stress  intensity  factor  ratio  :  Kt/K n  =  2230/0. 

Figure  54  :  Deformed  and  undeformed  meshes.  Inner  radius  1  mm,  outer  radius  10  mm. 
Stress  intensity  factor  ratio  :  Kj/Kji  =  1683/1252. 

Figure  55  :  Deformed  and  undeformed  meshes.  Inner  radius  1  mm,  outer  radius  10mm. 
Stress  intensity  factor  ratio  :  K|/Kj|  =  1098/1633. 

Figure  56  :  Deformed  and  undeformed  meshes.  Inner  radius  1  mm,  outer  radius  10mm. 
Stress  intensity  factor  ratio  :  K|/K||  =  0/2018. 

Figure  57  :  Strain  energy  density  along  a  circular  path  around  the  crack  tip  for 
selected  cases  of  mixed-mode  loading.  Radius  r=  1  mm  from  crack  tip. 

Figure  58  :  Full  logarithmic  representation  of  the  strain  energy  density  along  the 
iine  9=0°  for  selected  mixed-mode  cases  . 

Figure  59  :  Graphical  illustration  of  the  interval  halving  method. 

Figure  60  :  Angular  variation  of  the  stress  functions  crrr  ,cr 6Q  ,cre  and  crr0 

Stress  intensity  factor  ratio  :  K(/K||  =  2230/0. 

Figure  61  ■  Angular  variation  of  the  stress  functions  <jrr  ,<Tqq  ,cre  and  are 
Stress  intensity  factor  ratio  :  K|/K n  =  2222/175. 


Figure  62  :  Angular  variation  of  the  stress  functions  <j„  ,Cqq  ,cre  and  crrQ 
Stress  intensity  factor  ratio  :  K|/K||  =  2107/670. 

Figure  63  :  Angular  variation  of  the  stress  functions  crrr  ,aQQ  ,crQ  and  <7>g 
Stress  intensity  factor  ratio  :  K|/K||  =  1927/987. 

Figure  64  :  Angular  variation  of  the  stress  functions  a rr  ,< x00  .a e  and  crr0 
Stress  intensity  factor  ratio  :  K|/Kj|  =  1683/1252. 

Figure  65  :  Angular  variation  of  the  stress  functions  ,<?00  -O’e  anc*  ^>e. 

Stress  intensity  factor  ratio  :  Kj/Kj|  =  1405/1462. 

Figure  66  :  Angular  variation  of  the  stress  functions  crrr  ,cr 00  ,cre  and  c rr0 
Stress  intensity  factor  ratio  :  K(/K||  =  1098/1633. 

Figure  67  :  Angular  variation  of  the  stress  functions  cr„  ,cree  ,ae  and  ar6 
Stress  intensity  factor  ratio  :  K|/Kj|  =  772/1774. 

Figure  68  :  Angular  variation  of  the  stress  functions  a „  ,a ee  and  a:Q 
Stress  intensity  factor  ratio  :  K|/K|j  =  396/1903. 

Figure  69  :  Angular  variation  of  the  stress  functions  crrr  ,ae  and  arQ 
Stress  intensity  factor  ratio  :  K|/K|(  =  0/2018. 

Figure  70  :  Ratio  of  cree/crrg  along  the  line  9  =0°  for  selected  mixed-mode  cases. 

Figure  71  :  ln  versus  plastic  mixity  parameter  MP1  for  all  mixed-mode  cases 
considered.  Comparison  with  the  values  of  ln  obtained  by  Shih  [27]. 
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Table  1  :  Kj  -  K||  -values  according  to  the  fracture  criterion  by  Sih. 


Table  2  :  Shape  functions  of  the  eight  node  isoparametric  plane  strain  element. 


Table  3  :  Chemical  composition  and  material  data  of  the  stainless  steel  A304. 

Table  4  :  Comparison  of  J-integral  values  obtained  by  the  virtual  crack  extension 

method  and  the  direct  integration  method. 


Table  5  :  Powers  characterizing  singular  behavior  of  the  effective  von  Mises  stress 

along  the  line  0  =  0°  ahead  of  the  crack  tip. 


Table  6  :  Fracture  angle  0  and  corresponding  strain  energy  density  0.4  mm  from  the 

crack  tip  for  all  cases  of  mixed-mode  cases  considered. 

Table  7  :  Powers  of  the  singularity  of  the  strain  energy  density  in  the  vicinity 

of  the  crack  tip  along  the  line  0=0°  from  least  square  approximation. 


crack  length 
flow  vector 

nodal  displacement  vector 

slope  of  the  uniaxial  stress-strain  curve 

vector  of  body  forces 

strain-displacement  matrix 

element  mapping  matrix 

constant  term 

radius  of  the  HRR  field 

linear  elastic  constitutive  matrix 

elastic-plastic  constitutive  matrix 

Young's  modulus 

stress  functions  for  linear  elastic  material  behavior 
yield  criterion 
nodal  force  vector 

nodal  force  vector  caused  by  body  forces 
fiodal  force  vector  caused  by  surface  tractions 
function 

displacement  functions  for  linear  elastic  material  behavior 

global  energy  release  rate 

energy  release  rate  in  x  -  and  y  -direction 

constant 

J-integral  refering  to  crack  extentions  in  x  -  and  y  -direction 
J-integral  vector 

resultant  J-integral  value  of  the  J-integral  vector 
Jacobian  matrix 


ki .  k2 


stress  intensity  factor  refering  to  mode  I  and  mode  II 


K  |  =  k1  V  ji  stress  intensity  factor,  mode  I 

K||  =  k2 Vjt  stress  intensity  factor,  mode  II 

K  element  stiffness  matrix 

Ky  tangential  stiffness  matrix 

I  length 

n  strain  hardening  index 

Nj  shape  fanctions  for  element  nodes 

N  shape  function  matrix 

r  radius 

r0  inner  radius  of  the  HRR  field 

s  distance 

S  vector  of  surface  forces 

S  quadric  strain  energy  density  functional 

I  vector  of  surface  tractions 

u,  v  displacement  components  in  x  -  and  y  -direction 

U  displacement  vector 

u  displacement  function 

S  quadric  strain  energy  density  functional 

V  volume 

W  strain  energy  density  (general) 

We  external  work 

W(  internal  strain  energy  density 

x,  y,  z  cartesian  coordiantes 


material  constant 

vector  of  polynoiai  coefficients 

angle  of  crack  inclination 

exponent 


r 


e  ,  r) 


Subscripts 


boundary 

interpolation  polynomial 

vector  of  Uj  displacements  of  an  element 

strain 

vector  of  strain  components 

angle  (general) 

angle  of  crack  extention 

hardening  parameter 

shear  modulus 

residual  force  vector 

stress 

vector  of  stress  components 
stress  functions 
shear  stress 
Poisson’s  ratio 
intrinsic  coordinates 


min,  max 


Superscripts 


components  referring  to  cartesian  or  polar  coordinate  system 

critical 

elastic 

plastic 

maximum,  minimum 


transpose 


Historically,  conventional  stress  analysis  was  based  on  the  assumption  of  flawless 
material  behavior.  Since  the  existence  of  crack-like  flaws  cannot  be  precluded  in  any 
engineering  material,  fracture  theories  had  to  be  developed  which  account  for  local  stress 
concentrations. 

The  significance  of  intense  and  localized  concentration  of  stresses  around  sharp  notches 
was  first  emphasized  by  Inglis  [1].  He  realized  through  considerations  of  the  stress 
concentration  around  an  elliptical  hole  that  the  stress  becomes  infinitely  large  at  the  tip  of  a 
sharp  crack.  Based  on  the  ultimate  stress  concept,  this  would  indicate  that  a  cracked 
component  cannot  sustain  any  loading. 

Griffith  [2,3]  applied  energy  conservation  principles  to  the  problem  of  a  cracked  glass 
plate.  This  work  set  the  theoretical  foundation  for  the  field  of  continuum  fracture  mechanics. 
Irwin  [4]  and  Orowan  [5]  subsequently  modified  the  original  Griffith  theory  so  that  it  could 
be  applied  to  metals  by  adding  a  term  involving  the  plastic  energy  dissipation  rate  in  the 
plastically  deformed  region  near  the  crack  tip.  Due  to  difficulties  in  the  practical 
application  of  the  energy  balance  concept,  new  approaches  had  to  be  found  to  characterize  the 
material  behavior  under  the  influence  of  a  sharp  crack.  Irwin  [6]  was  able  to  utilize  the 
cracked  body  solutions  of  Westergaard  [7]  to  establish  a  relation  between  the  strain  energy 
release  rate  G,  (a  global  parameter)  and  the  stress  intensity  factor  K  (a  local  crack  tip 
parameter).  These  stress  intensity  factors  can  be  related  to  three  independent  local 
movements  as  shown  in  Figure  1.  These  are  categorized  as  : 

-  Mode  I,  or  opening  mode 

-  Mode  II,  or  sliding  mode 

-  Mode  III,  or  tearing  mode. 

Any  crack  deformation  in  the  case  of  linear  elastic  material  behavior  can  be  idealized  by  the 
appropriate  superposition  of  theses  cases.  Unlike  the  brittle  glass  considered  by  Griffith, 
most  metals  exhibit  the  phenomenon  of  ductility.  Crack  tips  are,  therefore,  engulfed  by 
plastic  yield  zones  with  finite  stresses. 

Early  attempts  to  model  the  plastic  deformation  surrounding  the  crack  tip  were  based 


upon  extensions  of  the  linear  elastic  fracture  mechanics  (LEFM).  Irwin  [81  broadened  the 
applicability  of  LEFM  by  introducing  a  modified  stress  intensity  factor  Kp.  At  the  same  time 
Wells  [91  established  the  crack  opening  displacement  (COD)  as  a  parameter  governing  crack 
extension  even  beyond  general  yielding.  Dugda'e  [10)  extended  the  COD  approach  and 
established  a  relation  between  a  plastic  zone  estimate  around  the  crack  tip  and  the  crack 
opening  displacement  in  thin  sheets. 

A  significant  contribution  to  the  field  of  elastic-plastic  fracture  mechanics  (EPFM)  was 
the  introduction  of  the  path  independent  J-integral.  This  integral  (originally  derived  by 
Eshelby  [ill  and  Cherepanov  [12))  was  introduced  into  the  field  of  fracture  mechanics  by 
Rice  [131.  Begley  and  Landes  [1 4]  showed  its  applicability  as  a  parameter  describing  the 
stress  concentration  at  the  crack  tip  and  suggested  the  use  of  a  critical  J-integral  value  J|c 
to  predict  the  onset  of  stable  crack  growth. 

Several  attempts  have  been  made  in  recent  years  to  arrive  at  a  more  general  definition  of 
the  J-integral  which  would  minimize  the  assumption  of  elastic  material  behavior  and  the 
absence  of  body  forces  while  still  retaining  its  desirable  features.  Some  of  the  proposed 
formulations  extend  the  definition  of  J  to  axisymmetric  three  dimensional  problems,  others 
consider  more  general  loading  conditions  (15-20). 

Hutchinson  [21]  and  Rice  and  Rosengren  [22]  independently  determined  the 
characteristic  singular  behavior  of  stresses  and  strains  inside  the  plastic  zone  (using 
deformation  theory  of  plasticity)  where  elastic  strains  are  negligible  compared  to  plastic 
strains.  This  zone  is  commonly  referred  to  as  the  HRR  singular  field  due  to  its  distinct 
singular  character  in  terms  of  stresses  and  strains.  In  their  analysis,  which  took  full 
advantage  of  the  path  independence  of  the  J-integral,  they  showed  that  stress,  strain  and 
displacement  components  can  be  related  to  dimensionless  functions.  These  functions  are  only 
dependent  on  the  hardening  characteristics  of  the  material  and  whether  the  material  is  in  a 
state  of  plane  stress  or  plane  strain.  Stress,  strain  and  displacement  components  inside  the 
HRR  field  are,  therefore,  determined  by  an  appropriate  stress,  strain  or  displacement 
function  and  a  singular  term  involving  the  J-integral  value  which  characterizes  the 
amplitudes  of  these  fields.  The  validity  of  expressing  crack-tip  quantities  in  terms  of  the 
HRR  singular  solution  has  been  shown  by  a  number  of  scientists  [23-26].  Shih  [24] 
established  (through  considerations  of  the  displacement  function  of  the  HRR-theory)  a 


relationship  between  the  J-integra!  and  the  COD-concept  and  proved  the  similarity  of  both 
concepts.  This  analysis  assumes  that  the  HRR  field  dominate  the  region  around  the  crack  tip 
having  a  size  of  at  least  ten  times  larger  the  crack-tip  opening  displacement. 

Shih  [27.28]  applied  the  HRR  singular  field  solution  to  the  case  of  a  body  under  mixed¬ 
mode  loading.  The  stress,  strain  and  displacement  functions  in  this  case  depend  on  the 
relative  composition  of  mode  I  and  mode  II  directly  ahead  of  the  crack.  The  J-integral,  in 
combination  with  a  parameter  sensitive  to  the  composition  of  mode  I  and  mode  II,  governs  the 
amplitude  of  the  singular  field. 

Though  many  cracks  in  structures  may  be  initially  under  mixed-mode  conditions,  most 
research  in  the  field  of  fracture  mechanics  is  focused  on  the  study  of  fracture  behavior 
under  pure  mode  I  conditions.  The  major  reason  for  this  is  the  general  observation  that  a 
crack  subjected  to  mixed-mode  loading  tends  to  grow  toward  a  mode  I  condition.  The  main 
interest  in  mixed-mode  fracture  mechanics,  therefore,  is  focused  on  determining  criteria 
which  predict  the  onset  of  crack  growth  and  the  angle  of  crack  extension  in  relation  to  the 
existing  crack.  In  contrast  to  pure  mode  I  where  the  critical  value  Jic  has  been  employed  to 
predict  the  onset  of  crack  growth,  this  quantity  is  no  longer  valid  for  mixed-mode  conditions 
since  the  crack  usually  does  not  extend  in  its  own  plane.  The  most  promising  concepts  of 
mixed-mode  fracture  criteria  are  therefore  based  on  energy  principles,  i.e.  the  maxima  or 
minima  of  either  the  total  strain  energy  density  or  of  its  components  [29-32].  Both  the 
concept  of  the  strain  energy  criterion,  introduced  by  Sih  [29]  and  the  T-criterion  suggested 
by  Theocaris  [33]  have  been  extended  for  use  in  the  elastic-plastic  regime.  Since  these 
criteria  are  of  local  nature,  they  depend  on  the  local  stress  and  strain  response  of  the 
material. 

In  the  analysis  of  cracked  bodies,  the  finite  element  method  has  become  the  major 
numerical  technique  for  the  solution  of  fracture  problems  (both  linear  and  nonlinear).  The 
theory  of  incremental  plasticity,  which  is  usually  incorporated  in  modern  finite  element 
programs,  relates  incren  -  its  of  stress  to  increments  of  strain.  The  formulation  of  the 
incremental  theory  of  plasticity  accounts  for  elastic  unloading  effects  and  has  been  very 
successful  in  simulating  ductile  material  behavior. 

In  the  present  investigation  local  crack-tip  field  quantities  were  examined  for  mixed¬ 
mode  loading  ranging  from  pure  mode  I  to  pure  mode  II.  A  ductile  material  was  modeled  and 
the  commercial  finite  element  package  ABAQUS  was  employed  to  perform  the  calculation.  All 


considerations  were  based  on  the  assumption  of  small  scale  yielding  which  requires  the 
plastic  zone  to  be  small  relative  to  other  dimensions  (i.e.  crack  length  or  specimen  size). 

Although  a  variety  of  mixed-mode  fracture  specimens  have  been  suggested  in  the  past 
[34,35],  there  is  still  no  universally  accepted  standard  mixed-mode  specimen.  In  this 
work,  therefore,  the  local  crack-tip  region  is  modeled  without  employing  a  mixed  mode 
fracture  specimen.  Displacements  on  the  boundary  are  calculated  by  assuming  elastic  stress 
intensity  factors  for  both  mode  I  and  mode  II  a  priori.  These  displacements  are  applied  to  the 
boundary  of  the  local  crack-tip  region. 

Ten  loading  combinations  which  span  the  range  from  pure  mode  I  to  pure  mode  II  were 
investigated.  The  J-integral  values  were  calculated  to  measure  the  strength  of  the  field 
singularities.  Stresses,  strains  and  both  total  and  elastic  strain  energies  were  examined 
with  respect  to  their  singular  behavior  and  angular  variation  within  the  plastic  zone. 
Details  of  plastic  zone  sizes  and  shapes  as  well  as  the  crack  blunting  under  varying  mixed 
mode  I  and  II  contributions  were  investigated.  The  stress  functions  for  all  investigated  cases 
were  determined  as  well  as  parameters  describing  the  amplitude  of  the  plastic  near-tip 
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The  stress  and  displacement  fields  around  a  crack  in  a  linear  elastic  material  have  been 
investigated  by  a  number  of  scientists.  Although  a  basic  solution  was  worked  out  by 
Muskhelishvili  [36],  Westergaard  [7],  Williams  [37],  Irwin  [6]  and  Sih  [38]  solved  the 
same  problem  using  different  approaches.  A  decisive  step  in  linear  elastic  fracture 
mechanics  (LEFM)  was  the  introduction  of  the  stress  intensity  factor  K  by  Irwin  [6].  His 
work  employed  using  Westergaard's  solution  for  the  near-tip  stress  field  of  a  cracked  body. 
If  the  elastic  solution  for  a  cracked  body  is  available,  the  stress  intensity  factors,  K|  and  K||, 
can  be  defined  as : 


K|  -  lim  crVY(r,0-O)(2jtr)i/2 

r—  0  ” 

(2.1) 

Kn  -  lim  crxy(r,0-O)(2Kr)i/2. 
r— 0  1 


Both  stress  and  displacement  fields  are  based  on  the  linear  theory  of  elasticity  and  may, 
therefore,  be  superimpose  1  The  stresses  and  displacements  under  combined  mode  I  and  mode 
II  can  be  written  as  for  the  coordinate  system  (given  in  Figure  2  )  : 
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where  E  is  Young's  modulus, 
u  is  Poisson's  ratio  and 

k  is  given  for  the  case  of  plane  strain  as  k  =  3-4u 


By  expressing  equation  (2.2)  in  the  form  : 


U|'  TZ 9''<e,t  TTva^  «i"(9) 
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several  characteristics  of  these  solutions  can  be  observed.  The  stress  intensity  factors 
depend  only  on  the  applied  loading  and  consequently  determine  the  intensity  of  the  local  field. 
The  remaining  terms  depend  on  the  spatial  coordinates  around  the  crack  tip  and  determine 
the  distribution  of  the  field.  These  subdivide  into  a  singular  1/Vr  component  and  an  angular 
component  expressed  by  the  geometric  functions  fjj*,  fy”,  gj1  and  g^. 

Higher  order  terms  of  the  actual  series  solution  have  been  neglected.  More  higher 
ordered  terms  need  to  be  included  if  the  field  had  to  match  outer  boundary  conditions.  Eftis  et 
al.  [39,40],  in  revisiting  the  stress  and  displacement  fields  of  the  one  parameter 
representation  given  in  equation  (2.2),  proved  the  inaccuracy  of  these  relations  for  the 
case  of  biaxial  loading.  This  stems  from  the  arbitrary  omission  of  the  second  term  of  the 
series  expansion  for  the  stress  components  which  contain  a  term  independent  of  the  distance 
from  the  crack  but  dependent  of  tne  angular  position  around  the  crack  tip.  Eftis  improved  the 
singular  solution  for  the  inclined  crack  under  biaxial  tensile  load  by  including  this  term 
which  finally  affects  only  the  x-component  of  the  stress  field.  Theocaris  et  al.  [41] 
developed  a  closed  form  solution  solution  of  stresses  and  displacements  of  a  slant  crack  under 
biaxial  tensile  loading  for  arbitrary  radius  (r)  away  from  the  crack  tip. 

Irwin  [42]  derived  the  relationship  between  the  stress  intensity  factors  and  the  energy 


release  rate  (G)  for  cracks  extending  in  their  original  plane.  For  pure  mode  I  and  pure  mode 
II  (under  plane  strain  conditions)  these  relations  are  given  as  : 


*1  2 
Gr  —  (1  -  u  ) 


<v  ~ 


(1  -  u2) 


(2.4) 


Energy  release  rates  can  be  added  for  a  crack  remaining  in  its  plane  according  to  : 

G  =  G|  +  G||  (2.5) 

For  a  linear  elastic  material  subjected  to  pure  mode  I  loading  conditions  the  energy 
available  to  create  a  unit  surface  is  G| .  The  critical  strength  parameter  governing  failure, 

therefore,  can  be  expressed  as  G|c-  This  parameter  can  be  related  to  the  critical  stress 
intensity  factor  KiC  by  equation  (2.4). 
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In  order  to  assure  a  comparison  of  field  quantities  for  arbitrary  ratios  of  stress  intensity 
factors,  a  criterion  had  to  be  applied  which  relates  Kt|-values  to  a  given  K| -value.  In 

contrast  to  a  pure  mode  I  Griffith  crack,  a  mixed  mode  crack  does  not  necessarily  extend  in 
its  own  plane.  Since  the  direction  of  crack  extent  is  not  known  a  priori  it  would  be  incorrect 
to  obtain  the  mixed  mode  energy  release  rate  by  adding  G|  and  G||.  In  contrast  to  pure  mode  l 
fracture  analysis  where  the  fracture  criterion  is  founded  on  a  given  K|c-value,  there  is  sill 
no  well  established  criterion  for  the  case  of  mixed-mode  fracture.  The  most  widely  used 
mixed  mode  fracture  criteria  are  [33,43-45]  : 

-  criterion  of  maximal  tangential  stress 

-  various  criteria  based  on  the  energy  release  rate 
various  criteria  based  on  the  energy  density 

-  criterion  according  to  Di Leonardo 

-  principal  strain  criterion 

-  J-integral  criterion 

-  modified  T-criterion  and 

-  various  empirical  criteria  based  on  experimental  K|-K||  failure  curves. 

The  energy  density  criterion  introduced  by  Sih  [29]  was  chosen  for  two  reasons: 

it  is  generally  in  good  agreement  with  experimental  results  [29,47]  and 
it  provides  a  concise  relationship  between  the  critical  energy  density  factor  and  the 
stress  intensity  factors  for  mode  I  and  II. 
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Sih  [29]  proposed  a  criterion  based  on  the  strain  energy  density  in  the  vicinity  of  a 
crack  tip.  For  an  elastic  material  the  strain  energy  density  is  given  in  its  general  form  as  : 


(~)  -  ^(<rx2+cry2+<722)-^((rx<7y+cryo-2+£rzo-x)+~(Txy2+rxy2+Ty22) 


where  E  is  Young's  modulus, 
u  is  Poisson's  ratio  and 
u  is  the  shear  modulus. 


Substituting  the  stress  components  from  the  asymptotic  linear  elastic  two  dimensional 
stress  solutions  given  in  equation  (2.2),  the  strain  energy  density  can  be  obtained  as: 


(•“)-  y- --p-(a-nki2  +  2a12kik2  +  a22k22) 


(3.2) 


The  coefficients  are  given  as: 


an  -  — —  [(  3-4u-cos0)  (1+cos0  )] 

1  6U 

ai2-  —  sin©  (cos0-(1-2u)] 

8  u 

a22  .  -1—  [4  (1-u  )  (l-cos0)  +  (1  +cos0  )  (3  cos0  -1)]  . 

1  6U 

It  can  be  seen  that  the  strain  energy  density  is  characterized  by  a  1/r  singular  term 
where  r  is  the  distance  from  the  crack.  The  quadric  term,  S,  in  equation  (3.2)  can  be 
considered  as  a  material  constant  [29]  and  varies  only  in  the  angle  0  around  the  crack  tip, 

t 

Figure  2. 

Determination  of  pairs  of  Kr  K(J  -  values  with  respect  to  an  assumed  maximum  value  of 


s.-V  ■" 


i 


Determination  of  pairs  of  K|-  K||  -  values  with  respect  to  an  assumed  maximum  value  of 
K|c  is  based  on  taking  advantage  of  two  hypothesis  formulated  by  Sih  [29J: 


1 .  A  crack  will  extend  in  the  direction  of  maximal  potential  energy. 

2.  The  critical  intensity  Scr  of  this  potential  field  governs  the  onset  of  crack 
propagation 


The  potential  energy  per  unit  volume,  P,  is  defined  as: 


P.-f 


Therefore,  P  assumes  a  maximum  if  the  following  relations  hold: 


£.0  ;  A<0 


(3.3) 


at  0  -  e0 


(3.4) 


The  formulation  of  the  stress  intensity  factors  kj  and  k2  for  a  crack  inclined  by  an  angle  0, 
and  of  length  2a  under  tensile  stress  <r  (  see  Figure  3  )  is  given  as: 


k1  -  <j  Va  sin2  0 


(3.5) 


k2  -  cr  Va  sin  3  cos  0  . 

Substituting  these  expressions  into  equation  (3.2)  yields  an  expression  for  S  : 

S  -  k1 2  (  a1 1  sin  0  +  2  a1 2  sin  0  cos  0  +  a^  cos20  )  sin  0  .  (3.6) 

Again,  according  to  the  first  hypothesis  and  equation  (3.3),  S  has  to  be  a  minimum  if  P 
shows  a  maximum.  Differentiating  equation  (3.6)  with  respect  to  0  and  setting  the  result  to 
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2  (1  -2u )  sin  (80  -  8)  -  2  sin  [2  (0O  -  8)]  -  sin  280  =  0 


(3.7) 


The  critical  values  of  K|  and  K||  lie  on  a  curve  in  the  K|-K||  plane  determined  from  equations 
(3.5)  and  (3.7). 

A  FORTRAN  program  has  been  written  to  perform  the  outlined  procedure  to  find  pairs  of 
K|-Km  values  for  a  given  pure  mode  stress  intensity  factor  K|.  The  material  data  are 

presented  in  chapter  6.2.  For  values  of  8  ranging  from  0  to  tt/2  equation  (3.7)  was  solved 
numerically  using  the  Newton-Raphson  Method.  Figure  4  shows  the  plot  of  KM  values  over  K( 
for  an  assumed  pure  mode  I  value  of  K|  -  2230  N//mm3/2  . 

Ten  pairs  of  stress  intensity  factors  which  span  the  range  from  pure  mode  I  to  pure  mode 
II  are  given  in  Table  1  and  will  be  referred  to  in  all  further  considerations.  These  pairs  of 
stress  intensity  factors  represent  points  on  the  K|-K|t  curve  which  are  the  endpoints  of 
equal  length  curve  segments. 


A  characteristic  ot  plastically  deformable  materials  is  that  a  particular  material  can 
undergo  different  histories  of  response  prior  to  the  body  reaching  its  end  state. 
Reversibility,  therefore,  cannot  be  guaranteed  after  load  removal.  The  final  strain  is  found 
to  be  dependent  on  the  history  of  loading,  in  addition  to  the  beginning  and  final  loading.  This 
means  that  plastic  behavior  is  a  path  function  and  requires  the  use  of  an  incremental  theory 
where  strains  are  integrated  over  the  strain  path  whenever  the  total  induced  strain  is  to  be 
determined.  A  common  approach  (employed  in  this  work)  is  the  incremental  theory  of 
plasticity. 

The  deformation  theory  of  plasticity  is  based  on  an  assumed  nonlinear  elastic  material 
response.  Plastic  strains  depend  only  on  the  current  state  of  stress  and  are  independent  of 
the  path  leading  to  this  state.  This  theory  (though  contrary  to  the  observed  nature  of  plastic 
behavior)  is  computationally  far  less  expensive  than  the  incremental  theory  and  is 
therefore  widely  used.  Figure  5  depicts  the  basic  difference  between  these  two  basic 
approaches  in  elastic-plastic  modeling. 

The  incremental  stress-strain  relationship  of  an  isotropic  strain  hardening  material 
can  be  derived  on  the  basis  of  the  following  relations  : 

a  yield  criterion 
a  yield  function 
a  flow  rule 

the  assumption  of  strain  rate  decomposition  and 
-  the  linear-elastic  constitutive  relation. 


i  )  Yield  Criterion  : 


Various  criteria  have  been  suggested  in  the  past  to  predict  the  onset  of  yielding  in  a 
material  subjected  to  loading.  The  Von  Mises  yield  criterion,  which  is  widely  used,  is 
based  on  the  assumption  that  yielding  occurs  when  the  second  invariant,  J2,  of  the 
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deviatoric- stress  -tensor  reaches  a  critical  value,  J2  =  2/3  cry2.  With  the  deviatoric 


stress  tensor  s.:  defined  as  : 


S . .  =  c  ■  •  *  T”  6  •  •  <r  1  l. 

ij  w  1 1  3  ij  °kk 


The  critical  value  of  J2  is  given  as  : 


(4.1) 


2  -.2 

s  .  s. .  =  —  cr 
ij  'J  3  y 


(4.2) 


The  effective  stress  may  be  written  as  : 


cr  =■  >TL[(  cr x-<r y)2+(cr y-<r z)2  +(crz-cry)2  +  6Txy2  +  6tx22  +  6rxz2  ]  .  (4.3) 


)  The  Yield  Function  : 

The  amount  of  hardening  of  an  isotropic  strain  hardening  material  can  be  expressed  by 
the  amount  of  plastic  work  which  is  : 


Wp  =  J  cr  ij  (d€jj)p  -  k 


(4.4) 


where  (d€jj)p  is  the  plastic  components  of  the  strain  differential  and 
k  is  the  strain  hardening  parameter. 


The  integral  is  path  (history)  dependent.  Like  the  equivalent  von  Mises  stress  given  by 
equation  (4.3),  the  equivalent  plastic  strain  increment  (de)pcan  also  be  obtained  from 
the  second  invariant  of  its  incremental  strain  tensor  (dejp  as  : 


(d€)0  ^4t  (d€ij)p  (d€jj’)p  }p 


(4.5) 


iii)  The  Flow  Rule  : 

The  flow  rule  governs  the  plastic  flow  after  yielding  and  can  be  derived  from  various 
yield  criteria  by  using  the  concept  of  a  plastic  potential  (f).  This  method  proposes  that 
the  incremental  strains  resulting  from  a  stress  tensor  a  jj  are  found  by  using  : 


3  f 

(deij/p  *  X  a 


(4.6) 


where  f  is  termed  the  yield  function  and 
x  is  a  constant. 


The  strain  history  and  its  current  magnitude  can  be  determined  by  the  yield  function  (f). 
If  the  J2  flow  rule  is  used  ,  (see  equation  (4.3))  the  yield  function  can  be  expressed  as  : 


f(o-ji)  =  (Sji  Sn)=  const 


(4.7) 


df=  —  dp- 
ecru  u 


(4.8) 


Three  cases  of  (df)  are  possible  : 


df  <  0  :  elastic  unloading  of  an  elastic  plastic  material  occurs 
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df  -  0  :  neutral  loading  of  an  elstic-plastic  material  and 

df  >  0  :  plastic  loading  of  an  elastic-plastic  material. 


For  the  case  where  (f)  is  taken  as  the  Mises  criterion,  given  in  equation  (4.2),  taking 
the  derivative  yields  : 


3  f  3  J2 

“Sii 

and  equation  (4.6)  simplifies  to  : 


(4.9) 


(d€jj)p  =  X  Sjj  (4.10) 

Fquations  (4.9)  and  (4.10)  are  referred  to  as  the  Prantl-Reuss  equations. 


i  v )  The  Strain  Rate  Decomposition  : 

During  an  stress  increment  the  resulting  strain  increments  can  be  split  into  their 
elastic  and  plastic  part  : 


-  dfie  +  d£p 


(4.11) 


v  )  Elasticity  : 

The  elastic  strains  can  be  related  to  the  deviatoric  and  hydrostatic  stress  components  by 
the  relations  : 


(€ij]  8 


J'j 


1  -2  u 
3  E 


(4.12) 
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Elastic  and  plastic  parts  can  be  added  and  the  complete  incremental  relationship  between 
stress  and  strain  for  elastic-plastic  deformation  can  be  obtained  as  : 


<d€  i  j)  e 


2  u 


1  -2 
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(4.13) 


It  can  be  shown  that  the  complete  incremental  elastic-plastic  stress-strain  relation  can 
be  written  as  derived  in  [48,49]  : 


dc=Depd£  (4.14) 

The  elastic-plastic  matrix  Dep  is  given  as  : 

Dsp-Q-  d  0  dT-r  ;do-£a  (4.15) 

a+  d  a 

where  A  is  the  local  slope  of  the  uniaxial  stress-strain  curve  and  can  be  gained 
from  the  stress-  strain  curve  of  the  given  material, 
a  is  the  flow  vector  which  is  a  partial  differential  of  the  yield  criterion 
with  respect  to  its  components  and  is  given  in  equation  (4.9). 

Q  <s  the  elasticity  matrix  having  the  form  for  the  case  of  plane  strain  : 


(1  -2u)  (  1  +  2u) 


1  -u 
u 

0  0 


u  0 
1  -u  0 
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Today  the  finite  element  method  is  firmly  established  as  a  standard  numerical  procedure 
for  the  solution  of  engineering  problems.  Its  versatility  is  based  on  the  following 
characteristics: 

irregular  geometries  can  be  modeled 

any  kind  of  boundary  condition  can  easily  be  formulated  and 

it  provides  sufficient  accuracy  for  many  engineering  purposes. 

Especially  in  the  field  of  fracture  mechanics,  the  finite  element  method  has  been  proven 
to  be  an  efficient  numerical  method  to  model  the  response  of  a  body  under  the  influence  of  a 
sharp  crack. 

The  basic  idea  behind  the  finite  element  method  is  to  divide  a  given  structure,  body  or 
region  into  a  number  of  elements.  The  elements  can  be  two  or  three  dimensional.  A  discrete 
number  of  nodes  situated  on  the  element  boundaries  connect  the  elements.  In  structural 
problems,  the  finite  element  method  solves  the  response  of  a  model  which  is  subjected  to  a 
given  load  by  determining  the  nodal  displacements.  A  set  of  interpolation  functions  (which 
are  referred  to  as  shape  functions)  uniquely  define  the  displacement  state  within  each 
element.  The  formulation  of  the  shape  function  depends  on: 

-  the  number  of  nodes  in  each  element  (order  of  element)  and 

the  number  of  independent  degrees  of  freedom  in  the  problem  considered. 

The  finite  element  formulation  for  a  continuum  can  be  obtained  by  taking  advantage  of  the 
formulation  of  the  principle  of  virtual  displacement  (which  is  a  special  case  of  the  principle 
of  virtual  work)  or  of  the  principle  of  minimum  potential  energy  (which  assumes  elastic 
body  behavior). 

In  the  following  very  brief  introduction  of  the  basic  finite  element  formulation,  only  the 
two  dimensional  case  will  be  considered. 


The  principle  of  virtual  displacement  states  that  equilibrium  is  obtained  if  the  work 
done  by  external  forces  We  on  an  arbitrary  virtual  displacement  field  equals  the  increase  in 

strain  energy  (W()  of  the  system  [48],  This  relation  can  be  expressed  in  a  variational  form: 

(5Wi»dWe.  (5.1) 

The  principle  of  virtual  work  can  be  formulated  as  the  volume  integrals  of  the  variations 
in  strain  energy  density  and  the  sum  of  variations  of  external  energies  resulting  from  body 
forces,  surface  tractions  and  point  loads.  Employing  matrix  notation,  the  variation  in 
internal  strain  energy  density  is  given  as: 

dWj  *  J  (<J£  )T SL  dV  (5.2) 

where  <5sJ  represents  the  variation  of  the  strain  vector  e  -  [  €x.  ey  Txy  ]T 
and 

<r  is  the  stress  vector  cr  -  [  crx  oy  rxy  ]T. 

The  variation  of  the  external  work  can  be  expressed  as: 

dWe  »  J(du)T&dv  +  J  (du)TSdr  +  2  (du  )T  fp  (5.3) 

where  <Ju  is  the  variation  of  the  displacement  vector  u  ■  [  ui-  U2  ]T  . 
ti  is  the  vector  of  body  forces  ja  =  [  b^ ,  b2  ]T  , 

S  is  the  vector  of  surface  tractions  &  *  [  St,  S2  ]T. 
r  is  the  boundary  where  surface  tractions  are  applied  and 
fp  are  nodal  forces. 
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The  finite  element  approach  is  based  on  the  assumption  that  displacements  within  an 
element  are  adequately  described  by  a  polynomial.  The  second  order  rectangle  has  eight  nodes 
and  its  interpolation  polynomial  approximation  to  the  displacement  field  is  assumed  in  the 


♦  ■  cxj  *  Q2  x  +  a3y  ♦  +  a5xy  ♦  a6y2  +  ayxy2  +agx2y  (5.4) 


In  order  to  assure  interelement  compatibility,  equation  (5.4)  must  be  complete  in 
terms  of  a  specific  power.  The  eight  constants,  can  be  evaluated  by  solving  a  set  of  eight 
simultaneous  equations  if  the  nodal  coordinate  are  inserted  into  equation  (5.4)  and  the 
displacements  equated  to  the  appropriate  nodal  displacements.  Performing  this  operation, 
equation  (5.4)  becomes  : 


Uj  -  [1.  *i  -  yi .  ,  Xjyi  ,  y;2  ,  xiyi2  ,  Xj2y|] 


(5.5) 


i>Ca 


(5.6) 


where  6  is  the  vector  of  u  displacements  of  an  element  £_  =  [  u1f  ug  ,  U3,.  ,,u8  ]T. 

Solving  equation  (5.6)  the  vector  of  constants  a  can  be  obtained  in  terms  of  nodal 
displacements  by  : 

a  =  (5.7) 


The  vector  &.  can  be  substituted  back  into  equation  (5.6)  and  : 
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is  obtained.  From  this  relation,  the  shape  function  matrix  M  can  be  obtained  as  : 
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u  = 


where  M  is  the  shape  function  matrix  N.  -  [N1t  N2 . N7,  N8)T 


(5.9) 


Similarly,  equation  (5.9)  holds  for  the  v  component  of  the  displacement  vector.  Equation 
(5.9)  can  easily  be  generalized  for  the  displacement  vector  u  by  writing  : 


u-  M  a 


(5.10) 


where  a  is  the  nodal  displacement  vector  a.  -  [ui,v1;  02,^2:  ■  ■  •  u7,v7.  u8,v8]T 
and 

N  is  the  shape  function  matrix  N.  -  [  i  N1 ,  [  N2 . 1  N7,  i  N8]T  . 

The  isoparametric  finite  element  formulation  has  proven  very  effective  in  structural 
analysis.  Isoparametric  elements  are  characterized  by  the  transformation  of  the  element 
geometry,  into  a  square  in  2  -  D  problems,  using  a  local  coordinate  system  defined  by  its 
5-  r\  coordinate  axes,  see  Figure  6  .  Axes  5  and  ri  pass  through  mid  points  of  opposite  sides, 
so  that  the  edges  are  defined  by  ?*±1  and  ri«±1.  If  the  shape  functions  used  to  describe  the 
geometry  and  displacements  of  an  element  are  the  same  then  this  element  is  called 
isoparametric.  The  shape  function  of  the  employed  isoparametric  eight  node  parabolic 
element  are  given  in  Table  2. 

The  displacement  components  of  any  point  within  the  element  are  defined  in  terms  of 
nodal  displacements.  In  equation  (5.2)  the  matrix  equation  for  the  strain  becomes  : 


s. 


(5.11) 


The  matrix  2.  is  defined  as  : 
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Assuming  linear  elastic  behavior,  the  stress-strain  relation  is  defined  through  the 
generalized  Hook's  Law  : 


sl  -  -  aaa 


(5.13) 


where  D.  is  the  linear  elasticity  matrix  given  in  equation  (4.15). 

In  the  isoparametric  finite  element  representation  the  shape  functions  Nj,  given  in  local 

coordinates  S,  n,  have  to  be  differenciated  with  respect  to  global  coordinates.  The  chain  rule 
must  be  applied  to  differentiate  : 


and 


3N  3N__3§_  9N  3n 
3x  *  3£  3x  +  3n  d* 


3n  3N__ae_  ±1 

3y  *  3?  dy  +  3ri  3y 


(5.14b) 


The  derivative  (3S/3x)  etc.  can  be  evaluated  from  the  inverse  of  the  Jacobian  matrix,  i*1. 
Using  the  Jacobian  matrix  the  volume  integral  (when  setting  dz-1  for  the  case  of  plane 
strain)  becomes  : 


dxdy  -  (detJ)d?  dn 


(5.15) 


Equation  (5.15)  can  be  substituted  into  equation  (5.2).  Employing  the  stress-strain 
relation  of  equation  (5.13)  and  employing  equations  (5.10)  and  (5.11)  the  finite  element 
formulation  (5.2) ,(5.3)  can  be  given  as  : 


{ if  fiT  (deti)  d?  dn  }  a  -  fb  +  fs  (5.16) 

where  fb  is  the  volume  integral  of  the  body  forces,  fb  =  11 NT  b  det  J.  d?  dn  and 
fs  is  the  integral  of  the  surface  tractions,  fs  =*  11 det  id?  dn. 

The  first  term  in  (5.16)  is  referred  to  as  the  element  stiffness  matrix  ft,  .  Equation 
(5.16)  can  be  numerically  integrated  using  the  Gauss-Legendre  quadrature  formula  where 
nine  integration  points  are  defined  for  the  isoparametric  plane-strain  rectangle,  see  Figure 
5.  The  integration  is  performed  in  the  ?,  n  space  where  the  coordinates  of  the  element  side 
range  from  -1  to  1. 


5.2  THE  FINITE  ELEMENT  FORMULATION  FOR  NONLINEAR  MATERIAL  BEHAVIOR 

The  described  finite  element  method  for  linear  elastic  material  behavior  can  be  extended 
to  materials  showing  nonlinear  behavior.  For  most  problems  in  material  plasticity  an 
incremental  algorithm  is  used.  It  is  based  on  the  incremental  theory  of  plasticity  where  the 
plastic  action  is  followed  as  it  develops,  and,  therefore,  accounts  for  the  path  dependence  of 
plasticity. 

The  initial  step  of  an  elastic-plastic  finite  element  calculation  assumes  linear  elastic 
behavior.  If  yielding  occurs  at  one  or  more  nodes  a  system  of  residual  forces  will  exist, 
such  that: 


-  ( I  +  J  NT  tl  dv  )  *  0 


(5.17) 


where  I  is  the  vector  of  applied  external  forces.  If  the  effective  stress  at  one  or  more  nodes 
exceeds  the  yield  stress,  the  material  stiffness  matrix  is  continually  varied.  Then 
increments  of  strains  are  related  to  increments  of  stresses  according  to  equation  (4.14)  : 

d  s.  =  Q^p  d£.  (5.18) 

where  is  the  elastic-plastic  matrix  given  by  equation  (4.15).  Equation  (5.18)  can  be 
substituted  into  (5.17)  and  a  relation  between  an  incremental  load  Au  and  the  increments  of 
the  residual  vector  AiL  (which  is  usually  not  zero)  is  obtained  as  : 

Ai  -  Kt  Au-(  A{  +  JnTAkdv  )*  0  (5.19) 


where  Ht  is  the  tangential  element  stiffness  matrix  in  the  elastic-plastic  range  and 
is  given  as  : 

K,  =  iaTc«padv . 

Equation  (5.19)  can  only  be  solved  iteratively  according  to  the  following  steps  : 

1 .  Employing  incremental  displacements  Ay.  in  each  iteration  step  r,  an  iterative 
correction  (dy)r  is  calculated  using  the  Newton  Raphson  Method  : 

(<?U)r  -  [  KTr  I'1  (5.20) 

2.  At  the  end  of  each  iteration  the  improved  displacement  estimate  is: 


Au.r+  =  Aur  +  (du)r 


(5.21) 


This  value  Au.r+1  is  substituted  in  (5.19)  to  evaluate  the  residual  force  vector  a^ 
which  is  used  in  (5.21)  to  calculate  the  correction  of  the  displacement. 

This  algorithm  is  repeated  until  the  maximum  of  the  residual  force  vector  is  smaller 
than  a  user  defined  number.  ABAQUS  uses  the  modified  Newton  Raphson  method  where  the 
stiffness  matrix  Kt  calculated  after  each  convergent  solution  instead  of  being  modified  after 
each  iteration.  This  results  in  a  significant  decrease  in  computing  time. 


5JL  GENERATION  OF  THE  FINITE  ELEMENT  MESH 

A  central  aspect  of  the  application  of  the  finite  element  method  is  the  generation  of  an 
appropriate  mesh.  The  quality  of  the  finite  element  mesh  affects: 

-  the  accuracy  of  the  solution 

the  amount  of  required  computing  time  and 
the  convenience  of  postprocessing  the  results. 

Today,  most  finite  element  meshes  are  generated  with  the  help  of  a  finite  element 
modeling  program.  For  the  present  investigation  the  mesh  of  the  mixed  mode  fracture  model 
was  created  on  an  IBM  5080  workstation  using  the  CAEDS  software  package  [49],  CAEDS  is  a 
computer  aided  design  tool  which  provides  the  ability  to  model  and  analyze  the  behavior  of 
mechanical  structures.  CAEDS  divides  this  task  into  three  consecutive  steps: 

1 .  Geometry  definition: 

The  model  geometry  is  defined  by  points  and  their  connecting  lines.  Subareas  have  to 

> 

be  defined  in  the  model  which  help  control  pattern  and  density  of  the  finite  element 
mesh  to  be  generated. 
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2.  Mesh  Generation: 

The  mesh  generation  accesses  the  geometry  of  the  model  through  defined  subareas. 
Concentration  and  configuration  of  the  finite  element  mesh  can  be  influenced  by 
defining  nodes  or  a  node  concentration  on  the  boundary  of  the  subarea.  Thus  the  finite 
element  mesh  for  every  subarea  is  generated  automatically  in  an  exactly  predictable 
manner. 

3.  Model  Checking: 

This  module  assures  the  correctness  of  the  created  mesh.  Internal  free  edges,  node 
and  element  coincidence  and  element  distortion  can  easily  be  detected  and  corrected. 
Furthermore,  the  bandwidth  and  the  profile  of  the  stored  matrix  can  be  optimized  and 
the  nodes  of  the  model  renumbered  accordingly  (which  shortens  the  computing  time 
significantly). 


v> 


For  the  present  investigation,  a  disk  shaped  fracture  model  with  a  sharp  crack  was 
modeled.  The  dimensions  of  the  specimen  were  mainly  determined  by  the  definition  of  small 
scale  yielding  which  limits  the  plastic  zone  size  to  approximately  20  percent  of  the 
specimen  size  [27].  Therefore,  for  the  given  sets  of  K(  and  K(|  values,  the  plastic  zone  size 
was  estimated  using  approximation  formulas.  Employing  the  plastic  zone  size  estimation 
formula  by  Irwin  [8]  yielded  a  maximum  plastic  zone  size  radius  of  rp*  3.75  mm  in  the 

case  of  pure  mode  l  (under  plane  strain  conditions).  The  yield  stress  contour  which  is  given 
in  Figure  7  for  a  selection  of  mixed  modes  indicates  in  the  case  of  pure  mode  II  loading  a 
plastic  zone  size  radius  of  approximately  rp»  24  mm.  The  specimen  radius  was  therefore 

chosen  to  be  100  mm.  The  crack  width  was  modeled  as  small  as  possible  to  produce  a  sharp 
crack. 

Conventional  elements  cannot  simulate  the  singularities  in  the  strain  fields  which  exist 
near  sharp  cracks  in  the  case  of  elastic  and  elastic-plastic  material  behavior.  Various 
authors,  therefore,  have  suggested  finite  elements  which  account  for  these  singularities 
without  using  large  numbers  of  elements. 

Henshell  and  Shaw  [50]  and  Barsoum  [51]  proposed  the  use  of  isoparametric  eight  node 
quadrilateral  elements  with  midsize  nodes  displaced  by  a  quarter  of  the  edge  length  towards 
the  crack  tip.  This  collapsed  quarter  point  element  produces  a  1/Vr  singularity  in  the 
elastic  strains.  Barsoum  [51]  proposed  that,  in  the  case  of  crack-tip  plasticity,  eight 
isoparametric  eight  node  elements  can  be  degenerated  into  triangular  shape  elements  by 
collapsing  the  element  at  their  crack-tip  nodes  without  shifting  the  midsize  node,  Figure  8. 
All  collapsed  nodes  at  the  crack  tip  remain  unconstrained  and  have  independent  degrees  of 
freedom.  It  has  been  shown  [52-54]  that  this  causes  three  effects  : 
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-  a  singularity  of  the  order  of  1/r  is  simulated  in  the  approximation  of  the  strain 
components.  This  coincides  very  well  with  perfectly  plastic  material  behavior  in  the 
near-tip  field 

-  the  ability  to  reproduce  large  strain  gradients  is  retained  and 

-  spurious  numerical  unloading  often  encountered  with  the  collapsed  quarterpoint 
element  is  eliminated. 

Comparison  with  analytical  results  performed  by  Shih  [52]  showed  that  this  element 
simulates  the  material  response  at  the  crack  tip  reasonably  well.  Barsoum  [51],  however 
has  shown  that  this  type  of  element  possesses  theoretically  unbounded  terms  in  the  stiffness 
matrix  but  which  are  usually  suppressed  by  the  smoothing  character  of  the  Gauss-Legendre 
quadrature. 

Figure  9  shows  the  finite  element  mesh  of  the  modeled  specimen  in  the  vicinity  of  the 
crack  tip.  A  fan  of  24  degenerated  elements  with  a  side  length  of  0.04  mm  defines  the  crack 
tip.  The  crack-tip  width  was  modeled  as  0.004  mm.  Adjacent  to  this  fan  is  an  intermediate 
zone  which  connects  the  crack-tip  fan  to  the  main  fan  consisting  of  23  circumferential 
layers  of  32  element  segments,  see  Figure  10.  Three  circumferential  element  layers  are 
needed  to  model  the  boundary  element  layer  of  16  elements,  (Figure  11),  where  the 
displacement  components  act  on  the  outer  nodes.  The  finite  element  model  of  the  specimen 
consists  of  1178  elements  and  3643  nodes.  In  order  to  investigate  the  accuracy  of  the  finite 
element  model,  the  elastic  stress  intensity  factors  were  calculated.  These  agreed  with  the 
input  values  to  within  0.1  percent. 


6.2  MATERIAL 


The  material  properties  of  the  stainless  steel  A304  [55]  were  used  as  the  material  data 
input  in  the  finite  element  calculation.  This  steel  finds  its  main  application  in  pressure 


containments  in  the  high  temperature  range  due  to  its  ability  to  sustain  high  plastic 
deformation  beyond  the  yield  stress.  Table  3  lists  the  material  data  of  the  steel  A304  and  its 
chemical  composition  given  by  Newman  [56].  The  stress-strain  behavior  of  the  employed 
material  was  modelled  using  the  Ramberg-Osgood  relation.  In  the  case  of  the  power 
hardening  simulation  of  the  material  response,  the  uniaxial  stress-strain  relation  is  given 
as : 


for  cr  <  cry 


for  cr  >  o-y 


(6.1) 


where  cr  is  the  uniaxial  tensile  stress, 

€y  is  the  yield  strain, 
cry  is  the  yield  stress, 

a  is  a  material  constant  which  is  given  as  0.75  for  the  steel  A304  and 
n  is  the  hardening  index  which  is  given  as  6  for  the  present  material. 


Figure  12  shows  the  modeled  stress-strain  behavior  of  the  employed  steel  A304. 


6.3  USING  THE  FINITE  ELEMENT  SOLVER  'ABACUS' 

The  commercially  available  finite  element  package  ABAQUS  [57]  was  used  for  model 
solution.  Nodal  coordinates  and  element  connectivity  generated  by  CAEDS  can  be  accessed 
through  a  universal  file.  A  FORTRAN  program  has  been  written  to  both  reformat  the 
universal  file  for  the  correct  ABAQUS  input  and  to  correct  the  node  numbering  direction 
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since  CAEDS  does  not  employ  a  consistent  node  numbering  direction  within  an  element.  The 
applied  boundary  displacements  for  all  cases  of  mixed  modes  considered  were  calculated 
according  to  equation  (2.2).  The  element  type  CPE8  (eight  node  parabolic  plane  strain 
isoparametric  element)  was  employed  in  this  model. 

Based  on  the  modeled  stress  strain  curve  of  the  employed  material,  discrete  values  of 
stresses  and  plastic  strains  had  to  be  specified  in  the  ABAQUS  input.  Element  sets  for  both 
data  output  and  graphic  display  of  deformed  meshes  and  contours  of  specified  variables  were 
defined. 

For  all  cases  considered,  the  model  response  for  the  elastic-plastic  material  behavior 
was  calculated  using  small  strain  theory.  ABAQUS  generates  increment  sizes  automatically 
and  assumes  a  maximum  number  of  six  iterations  per  increment.  This  usually  assures  good 
convergence  at  relatively  short  computing  time. 

While  good  convergence  was  obtained  for  cases  of  mixed  modes  with  either  predominant 
mode  I  and  mode  II  contributions,  the  following  four  cases  of  mixed  modes  had  to  be 
subdivided  into  separate  steps  of  increasing  pairs  of  K|  -  K|(  values  to  obtain  a  convergent 
solution  : 


-  K|  /  K||  - 

-  K|/K„- 

-  K  |  /  K„  ■ 

-  K,/K„- 


2107/670 

1927/987 

1683/1252 

1405/1462 


Four  increasing  pairs  of  combinations  of  K|  and  Kjj  values  (of  equal  ratio)  were  assumed 

which  resulted  in  good  convergence  for  each  step.  J-integral  values,  however,  cannot  be 
obtained  by  ABAQUS  if  a  calculation  is  subdivided  into  separate  steps  [57]. 

In  order  to  permit  an  efficient  postprocessing  of  ABAQUS  data,  a  FORTRAN  program  has 
been  written  which  reads  any  variable  from  the  data  output  file  : 

-  along  a  line  having  a  specified  angle  to  the  x  -  axis 
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The  J-integral,  which  was  originally  established  by  Eshelby  [11),  was  introduced  by 
Bice  [13]  into  the  field  of  fracture  mechanics.  Rice  showed  that  the  potential  energy 
release  rate  for  a  two  dimensional  crack  extending  in  its  plane  in  a  homogeneous  linear  or 
non-linear  elastic  material  was  equal  to  a  path  independent  integral.  Its  definition  is  given 
in  cartesian  coordinates  as,  (see  Figure  13)  : 

3u  j 

J  -  J  Wdx2  -  Tj  ds  (7.1) 

r  1 

where  W  is  the  energy  density,  defined  as  W  =  J  cr,j  d€jj 
T  is  an  arbitrary  path  around  the  crack  , 

T j  is  the  traction  vector  defined  according  to  the  outward  normal  n  along  r 
and  is  given  as  Tj  -  c*jj  nf . 

The  J-integral  is  well  established  as  a  parameter  which  describes  the  magnitude  of 
near-tip  stress  and  strain  fields.  Knowles  and  Sternberg  [58)  subsequently  generalized 
the  J-integral  to  be  a  vector,  J^ ,  corresponding  to  the  potential  energy  release  rate  in  any 

coordinate  direction  of  the  crack  extension.  Jk  is  defined  as: 


where  nk  denotes  the  unit  outward  normal  to  r,  lying  in  the  same  plane  of  the  crack. 

For  the  two  dimensional  combined  mode  I  and  mode  II  fracture,  only  the  Jr  and  J2- 

integral  definitions  need  to  be  considered.  Both  integrals  have  the  following  important 
properties: 

i  )  Path  Independence  : 

A  proof  of  the  path  independence  of  the  J-integral  can  be  found  in  [59], 
Finite  element  investigations  of  J-integral  values  obtained  for  paths  very  close 
to  the  crack  tip,  however,  (assuming  elastic-plastic  material  behavior)  showed 
significant  path  dependence  where  J-integral  values  approach  zero  very  rapidly. 
McMeeking  [60,61]  investigated  this  behavior  systematically  and  showed  that  it  can 
be  related  to  the  large  deformations  around  the  blunting  crack.  The  J-integral  values 
calculated  along  paths  more  than  5  to  10  times  the  crack  opening  away  from  the  crack 
tip  can  be  considered  as  path  independent.  Their  role  as  parameter  characterizing  the 
crack-tip  field  quantities  is  probably  retained. 

i  i )  Compatibility  with  Linear  Elastic  Fracture  Mechanics  : 

For  linear  elastic  behavior  J1  and  J2  are  equivalent  to  the  energy  release  rate  G  in  x1 
and  x2  direction,  respectively.  Hellen  et  al.  [62]  and  Blackburn  [63]  related  Jr  and 
J2-integral  values  to  stress  intensity  factors  for  a  two  dimensional  crack.  The 
relations  are  given  for  the  case  of  plane  strain  as  : 
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iii) Application  in  Elastic  Plastic  Fracture  Mechanics: 

The  J-integral  value  for  two  dimensional  crack  problems  has  been  used  by  many 
authors  to  predict  the  onset  of  crack  growth  initiation  in  cracked  bodies  both  for 
linear  elastic  and  elastic  plastic  material  behavior.  This  concept  was  introduced  by 
Begley  and  Landes  [14],  [64].  A  critical  value  of  J|C  in  pure  mode  I  for  plane  strain 

conditions  can  be  determined  by  a  standard  test  method  if  the  conditions  of  quasistatic 
loading,  negligible  body  forces,  monotonic  loading  and  stationary  crack  are  met. 
Kishimoto  et  al.  [IS],  in  their  interpretation  of  J1  an  J2  as  vector  quantities 
of  the  strain  energy  release  rate  for  crack  extension  in  the  two  dimensional  case, 
defined  a  resultant  vector  Jres.  Its  magnitude  is  given  as  : 

Jres  -  ^  Ji2  +  J22  (7.4) 


They  proposed  that  failure  in  mixed  mode  occurs  if  the  resultant  J-integral  equals  the 
critical  energy  release  rate: 

Jres^Gcr  (7.5) 


According  to  Bakker  [65]  this  criterion  has  found  little  experimental  verification  for 
the  case  of  mixed  mode  fracture. 
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Various  procedures  have  been  developed  in  the  past  to  calculate  the  J-integral  value.  A 
survey  of  different  methods  can  be  found  in  [66].  In  the  present  investigation  the  virtual 
crack  extension  method  and  the  direct  integration  method  were  employed. 


i  )  The  Virtual  Crack  Extension  Method  : 


This  method  originally  described  by  Parks  [67]  as  the  stiffness  derivative  method  is 
an  implemented  feature  in  most  commercial  finite  element  programs  like  ABAQUS.  In 
this  method,  the  potential  energy  release  rate  is  evaluated  directly  in  a  single  finite 
element  analysis  by  advancing  the  crack  tip  or  crack  front  by  a  small  amount  .  This 
small  advance  changes  the  stiffness  of  some  of  the  elements  in  the  mesh  and  the  change 
in  potential  energy  can  be  calculated.  Even  for  coarse  meshes  this  method  yields  very 
accurate  results  [53].  According  to  Nagtegaal  [66],  the  method  of  virtual  crack 
extension  is  particularly  accurate  if  collapsed  elements  at  the  crack  tip  are  used. 


i  i )  The  Direct  Integration  Method  : 

A  direct  way  of  calculating  both  J^and  J2 -integral  values  is  the  numerical 
integration  of  equation  (7.1)  using  discrete  data  points  (e.g.  from  a  finite  element 
analysis).  The  circumferentially  arranged  elements  in  the  present  mesh  suggested  the 
use  of  circular  paths  around  the  tip.  The  definition  of  Ji  and  J2  had  to  be  expressed  in 
polar  coordinates.  The  transformation  of  x  and  y  into  polar  coordinates  is  given  as: 

f 

x  =  rcos  9,  y  =  r  sin  0  (7.6) 
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The  path  increment  ds  can  be  expressed  as  : 


ds  =  rde  (7.7) 

The  displacement  derivatives  with  respect  to  x,  and  x2  assume  the  following  form, 
respectively  : 
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The  traction  vector  can  be  expressed  along  a  circular  path  as  : 


(7.8) 
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J1  and  J2,  calculated  along  a  circular  path  can  now  be  given  as  : 


Ji 


dUj  1  dUj 
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(7.11) 


?  d  U  j  1  d  U  j 

J  2  »  r  J  ( W  sm  e  -  Tj  (  sin  e  -r—  +  —  oose  - — )]  de 

.  *  3  r  r  eg 


Here,  For  the  case  of  elastic-plastic  material  behavior,  the  strain  energy  density  is 
the  sum  of  its  elastic  and  plastic  (dissipative)  part : 


w  -  We  +  WP| 


(7.10) 


A  FORTRAN  program  has  been  written  to  perform  the  numerical  calculation.  Stress 
components,  elastic  energy  density  as  well  as  plastic  dissipation  and  displacements  were 
read  from  the  ABAQUS  output  file  along  a  defined  path  by  a  data  post  processing  program. 
All  quantities  related  to  an  integration  point,  (e.g.  stress  components,  linear  energy 
density  and  plastic  dissipation)  were  interpolated  linearly  for  the  adjacent  element  node. 
The  derivatives  of  the  displacement  components  were  computed  using  a  sixth  order  finite 
difference  formula.  The  integration  along  a  path,  defined  by  85  equally  spaced  element 
nodes,  was  performed  numerically  using  the  Simpson's  second  order  integration  rule. 

For  all  mixed  mode  cases  considered  both  Ji  -  and  J2-integral  values  were  calculated 
for  radii  ranging  from  0.6  to  8.3  mm  from  the  crack  tip. 


7.1.  RESULTS  QF  THE  J1  AND'Jg  INTEGRAL  CALCULATIONS 

Figure  14  depicts  the  variations  of  the  J^-  and  ^-integral  values  for  all  mixed  modes 
considered  over  the  radial  distance  r  from  the  crack  tip.  J1 -integral  values  range  between 
9.8  N/mm  for  the  case  of  pure  mode  II  to  21.2  N/mm  for  the  case  of  mixed  modes  given  as 
K|/K||-  2222/172.  J2  values  are  generally  negative  and  range  between  J2  -  0  for  the 
cases  of  both  pure  modes  I  and  II  and  -14.9  for  the  case  of  K|/K||  -1405/1462. 

Good  path  independence  was  obtained  for  all  cases  considered.  For  increasing  mode  II 
contribution,  however,  greater  variations  in  the  J-j -integral  values  can  be  observed.  The 
maximum  deviation  reaches  9.3  percent  for  pure  mode  II  loading  in  comparison  to  a 


36 


r 


f 


<5 


deviation  of  1.12  percent  in  the  case  of  pure  mode  I. 

J2-integral  values  are  generally  less  accurate  and  therefore  show  more  path 
dependence  in  their  results.  Hellen  et  al  [62]  pointed  out  that  due  to  higher  overall 
displacement  gradients  in  the  x2  -direction,  the  associated  numerical  errors  are  large. 
This  can  be  seen  from  the  deviations  of  J2  -  integral  values  which  range  between  9  percent 
(for  the  case  of  K|/K M  »  2222/175)  and  17  percent  (for  the  case  of  K(/ K u  » 
1405/1463). 

Both  Jr  and  J2-  integral  values  show  little  variation  for  outer  paths  between  4.1  and 
8.6  mm  and,  therefore,  these  values  may  be  viewed  as  more  accurate. 

Table  5  lists  the  Ji -integral  values  for  the  outermost  path  from  both  the  direct 
integration  and  the  virtual  crack  extension  performed  by  ABAQUS.  As  pointed  out  in 
chapter  6.2,  the  Jt -integral  calculations  of  four  cases  of  mixed  modes  could  not  be 
calculated  by  ABAQUS.  Deviations  of  J-| -integral  values  gained  from  both  methods  are 
within  5.7  percent. 

In  all  further  investigations  the  integral  values  given  in  Table  4  will  be  used. 


The  dominant  singularity  solution  for  a  cracked  plate  in  a  power  law  hardening  material 
has  been  given  independently  by  Hutchinson  [21]  and  Rice  and  Rosengren  [22]  for  both  mode 
I  and  mode  II  stress  distributions.  The  solutions  which  are  known  as  the  HRR  singular  field 
were  generalized  by  the  solution  for  the  mixed-mode  stress  distributions  presented  by  Shih 
[27].  For  the  small-scale  yielding  case,  the  region  around  the  crack  tip  can  be  divided, 
(according  to  the  nature  of  singular  material  behavior),  into  three  distinct  areas  [59]  : 

-  the  far  tip  field 

-  the  near  tip  fieid  and 

-  the  intermediate  zone. 

Figure  15  identifies  these  areas. 


8.1  THE  FAR  TIP  FIELD 


At  distances  large  compared  to  the  plastic  zone  size  the  stress  and  strain  distribution  is 
dominated  by  the  1/Vr  singularity  from  the  linear  elastic  solution  for  the  stress  and 
displacement  fields.  A  measure  of  the  strength  of  the  singularity  is  the  path  independent 
J-integral  which  can  be  related  to  the  stress  intensity  factors  K(  and  K|(  according  to 

equation  (2.2).  A  convenient  definition  which  characterizes  the  relative  strength  of  Kt  and 
K ||  in  the  far  tip  field  was  introduced  by  Shih  [27]  as  : 
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(8.1) 


M®  is  referred  to  as  the  far  tip  field  mixity  parameter  which  ranges  from  0  to  1  with  M®»1 


for  pure  mode  I  and  M®  »  0  for  pure  mode  II. 


For  a  strain  hardening  material  which  can  be  described  by  a  power  law,  i.e.  the 
Ramberg-  Osgood  relation,  the  stress-strain  relation  based  on  the  deformation  theory  of 
plasticity  is  given  as  : 


1 + u  „  .  1+  2  u  .  3  /<X\ 
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n  -  1  e  . 


(8.2) 


where 


s  the  yield  stress, 


s  the  deviatoric  stress  tensor  given  in  equation  (4.1), 


s  the  effective  stress  given  in  equation  (4.3), 


s  the  yield  strain, 


s  Poisson's  ratio, 


s  Young's  modulus, 


s  a  material  constant  and 


s  the  hardening  coefficient. 


Large  plastic  strains  can  be  expected  in  the  near  field  so  that  (with  negligible  elastic 
strains)  equation  (8.2)  becomes  : 


n  -  i  e  •• 
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(8.3) 


It  can  be  assumed  that  the  only  singularity  contained  in  this  region  is  associated  with  the 
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crack  tip.  For  a  circular  path  of  radius  r,  where  r0  <  r  <  D,  enclosing  the  crack  tip  (see 
Figure  15  ),  the  J-integral  (according  to  equation  (7.11)  )  remains  path  independent. 

To  ensure  the  path  independence  of  the  J-integral  value  the  integrand  must  exhibit  a  1/r 
singularity.  Since  the  integrand  is  essentially  a  product  of  stress-  and  strain-like 
components,  this  product  leads  to  a  function  (f)  which  is  only  dependent  on  e  multiplied  by  a 
1/r  term  (assuming  the  material  behavior  satisfies  equation  (8.3))  : 


Hutchinson  [21]  has  shown  this  to  be  the  case  for  power  law  hardening  materials  if  the 
stresses  and  strains  are  given  in  polar  coordinates.  For  a  power  hardening  law  satisfying 
equation  (8.3),  equation  (8.4)  implies  that  the  following  relations  hold  : 
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where  cr  jj,  €jj  and  Uj  are  stress,  strain  and  displacement  functions  in  polar 

representation  where  i  and  j  are  radial  and  angular  components  and 
C  is  a  material  dependent  constant  term. 

Different  from  the  asymptotic  s.lution  for  linear-elastic  material  behavior,  therefore,  the 
singular  fields  in  the  elastic-plastic  range  are  dependent  on  the  hardening  characteristic  of 
the  material. 

Rice  and  Rosengren  [22]  and  Hutchinson  [21]  solved  equation  (8.5)  for  the  stress, 
strain  and  displacement  functions  by  introducing  an  Airy  stress  function.  A  partial 
differential  equation  governing  the  stress  function  can  be  derived  from  the  compatibility 
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equation  which  can  be  reduced  to  a  fourth  order  nonlinear  differential  equation  and  solved  by 
a  higher  order  finite  difference  scheme.  A  more  detailed  discussion  of  the  procedure  can  be 
found  in  [59]. 

The  constant  term  C  can  be  determined  by  taking  advantage  of  the  path  independence  of 
the  J-integral.  Substitution  of  equations  (8.5)  into  the  definition  of  the  J-integral 
(equation  (7.1))  leads  to  the  determination  of  the  constant  C  as  : 
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where  a  is  a  material  constant  , 

€y  is  the  yield  strain, 

cry  is  the  yield  stress, 

n  is  the  hardening  exponent  and 
ln  is  a  constant  given  by  : 
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These  equations  (which  were  originally  formulated  for  pure  mode  I)  can  be  extended  for 
the  case  of  mixed  modes.  Due  to  the  path  independence  of  the  J i  -integral  (regardless  of 
mixed-mode  contributions),  all  near  tip  field  quantities  remain  under  the  control  of  the 
J ! -integral  value. 

In  the  same  manner  as  the  elastic  mixity  parameter  (M®  in  equation  (8.1)  )  the  plastic 
mixity  factor  Mp  identifies  the  relative  composition  of  mode  I  and  mode  II  directly  ahead  of 
the  tip  according  to  : 
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From  equations  (8.5  to  8.7)  the  formulas  defining  stress,  strains  and  displacements  in 
the  near  field  of  a  crack  under  mixed-mode  conditions  can  now  be  given  as  : 
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The  stress  and  displacement  fields  are  therefore  characterized  by  a  1/r(1/(n^1 ) 
singularity  whereas  the  strain  field  assumes  a  f/r(n/(n+1)  singularity.  In  reality  such 
large  stress  components  cannot  exist  since  geometry  changes  modify  several  aspects  of  the 
tip  field  and  therefore  limit  the  stress  concentration  at  the  tip  (as  indicated  in  the  blunting 
analysis  of  McMeeking  [60,61]). 


8.3  THE  INTERMEDIATE  ZONE 

Combination  of  the  HRR  field  and  the  far  field  characterizes  the  stress  and  strain 
distribution  and  magnitude  of  the  intermediate  zone.  Whereas  its  outer  border  is  defined  as 
the  transition  from  the  elastic  to  the  plastic  zone,  its  border  to  the  HRR  field  can  not  be 
distinguished  clearly.  In  general  it  can  be  assumed  that  the  powers  characterizing  the  field 
singularity  of  the  intermediate  zone  show  a  smooth  transition  into  the  characteristic 
powers  of  the  HRR  field.  No  analytical  solution  has  been  found  yet  to  connect  near  tip  field 
and  far  tip  field  quantities  [27,59]. 
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Stresses  and  the  strain  energy  densities  around  the  crack  tip  of  the  investigated  plane 
strain  specimen  were  studied  with  respect  to  their  singular  behavior  under  varying  mixed¬ 
mode  conditions.  All  data  were  taken  from  the  finite  element  analysis  output. 


Figure  1 6  shows  the  plot  of  the  effective  stress  versus  the  distance  ahead  of  the  crack 
tip  along  the  line  0*0°  for  all  mixed-modes  cases  considered.  Four  important  features 
characterize  the  material  response  under  the  influence  of  mixed  modes  : 

i )  The  effective  stress  increases  sharply  for  higher  mode  II  contributions  under 
comparable  loads.  This  tendency  is  also  observable  if  the  effective  stresses  around  the 
crack  tip  is  considered.  Figure  17  shows  the  effective  von  Mises  stress  along  a  circular 
path  having  a  radius  of  r  »  0.4  mm  away  from  the  crack  tip. 

i  i )  For  three  cases  of  low  mode  II  contributions  there  is  a  distinct  transition  zone  between 
the  elastic  and  the  elastic-plastic  zones.  Both  zones  are  separated  by  the  yield  stress  of 
cTy  ■  265  MPa. 

i  i  i )  The  elastic-plastic  zone  shows  in  its  outer  region  an  increasing  influence  of  the 
intermediate  zone.  Plastic  strains  which  are  of  the  order  of  their  elastic  counterparts 
result  in  a  combination  of  the  HRR  field  and  the  far  field. 

iv)  At  distances  very  close  to  the  crack,  i.e.  less  than  0.5  mm.  the  singularity  governing 
the  effective  stress  behavior  is  weaker  than  the  singularity  characterizing  the  elastic 
material  behavior  and  can  therefore  be  attributed  to  the  HRR  field.  In  Figure  18.  using 
full  logarithmic  axes,  the  nature  of  occurring  singularities  are  shown.  For  distances 
close  to  the  crack  tip,  the  equivalent  stresses  of  all  mixed  modes  are  distinguished  by 
parallel  lines  before  the  smooth  transition  into  the  intermediate  zone  which  show 


particularly  weak  singular  behavior.  For  four  cases  of  strong  mode  I  contribution  the 
elastic  singularity  is  sharply  separated  from  the  elastic-plastic  zone. 

Figure  19  shows  the  y  -component  of  the  stress  tensor  versus  the  distance  from  the 
crack  along  a  line  9  «  0°  for  all  mixed-modes  cases  considered.  Increasing  mode  II 
contribution  results  in  a  steady  decrease  in  the  values  of  the  stress  component  in 
y-direction  which  approach  zero  for  the  case  of  pure  mode  II.  For  distances  less  than 
approximately  0.2  mm  away  from  the  crack  tip  significant  scatter  in  the  y-component  of 
the  stresses  for  mixed  modes  of  high  mode  II  contribution  can  be  observed.  It  is  assumed  that 
this  stems  from  the  influence  of  the  badly  distorted  crack  tip  which  is  exposed  to  an 
increased  rotation  as  mode  II  contributions  grow.  Since  stresses  are  a  second  order  quantity, 
(that  means  they  are  calculated  from  displacement  derivatives  obtained  from  the  finite 
element  solution)  this  effect  may  be  amplified. 

The  powers  of  the  singularities  can  easily  be  determined  if  the  stress  is  assumed  to 
follow  the  form  : 

<r9~Crio-e  (9.1) 

In  full  logarithmic  notation  of  equation  (9.1)  the  parameter  7  indicates  the  slope  according  to  : 

fog  (o-e)  =.  7  log  (r)  +  log  (C  o-e)  (9.2) 

The  extraction  of  the  exponent  7  has  been  performed  using  a  least  square  approximation  for 
all  mixed-modes  cases  considered.  Table  5  lists  these  powers  for  distances  from  the  crack 
which  contain  the  characteristic  singularity.  Compared  to  a  predicted  value  of  the  power  of 
7—1/7  for  the  employed  material,  it  can  be  seen  that  for  increasing  mode  II  contributions 
this  value  is  approached.  It  reaches  the  predicted  value  of  7  ■  -0.1428  almost  exactly  in  the 
case  of  pure  mode  II.  For  overwhelming  contributions  of  mode  I  the  HRR  solution 
characterizing  power  was  not  reached  even  for  closest  distances  to  the  crack  tip.  This  is 
mainly  caused  by  the  distinct  crack-tip  blunting  which  results  in  a  decrease  in  the  effective 
stress  and  weakens  their  singular  behavior.  The  elastic  singularity  which  could  be 
determined  for  four  cases  showed  excellent  agreement  with  the  predicted  value  of  7  ■  -0.5  . 


Contour  plots  of  the  effective  von  Mises  stress  in  the  vicinity  of  the  crack  tip  allow  a 
good  qualitative  assessment  of  the  body  response  under  the  influence  of  a  crack  under  mixed¬ 
mode  loading. 

Figures  20  to  26  display  the  contours  of  the  effective  stress  in  a  zone  of  1  mm  radius 
around  the  crack  tip  for  a  selection  of  mixed-modes  cases  (which  are  identified  on  each 
plot).  The  typical  symmetric  butterfly  shape  of  the  stress  contours  in  the  case  of  mode  I 
incline  and  assume  asymmetric  shapes  through  stages  of  mixed  modes  with  increasing  mode 
II  contributions  until  the  contours  show  the  typical  compact  and  symmetric  mode  II  pattern. 

The  increasing  concentration  of  the  von  Mises  stresses  around  the  crack  tip  is 
significant  for  higher  mode  II  contributions  While  for  overwhelming  mode  I  contributions 
elastic  regions  can  be  still  observed,  the  zone  considered  is  entirely  plastic  in  the  range  of 
higher  mode  II  values. 

The  shapes  and  magnitudes  of  the  effective  stresses  and  their  variations  between  1  mm 
and  10  mm  radii  for  mixed  modes  and  2  mm  to  20  mm  radii  for  pure  mode  II  around  the 
crack  tip  are  shown  in  Figures  27  to  33.  Here  the  influence  of  mode  II  contributions  results 
in  higher  effective  stresses  and,  consequently,  in  larger  plastic  zone  sizes. 

The  outermost  contours  of  the  equivalent  plastic  strains  represent  a  good  measures  of 
the  plastic  zone  size.  Figures  34  to  37  are  contour  plots  of  effective  plastic  strains  for 
selected  mixed-mode  cases  in  a  circular  region  of  1  mm  radius  around  the  crack  tip.  The 
increasing  gradient  of  plastic  strains  around  the  crack  tip  is  evident.  Figures  38  to  43  show 
the  effective  plastic  strains  in  the  region  of  1  to  10  mm  around  the  crack  tip  for  ail  mixed 
modes  and  the  region  2  to  20  mm  for  the  case  of  pure  mode  II. 

Figures  44  to  48  show  the  deformed  mesh  for  a  selection  of  mixed-modes  cases  in  a 
region  of  1  mm  radius  around  the  tip.  The  displacements  are  magnified  by  a  factor  of  two.  In 
the  case  of  pure  mode  I,  a  parabolic  shaped  crack  blunting  can  be  observed.  In  chapter  6.1  it 
was  pointed  out  that  the  employed  crack  tip  elements  can  only  approximate  the  strain 
assymptote.  For  increasing  influence  of  mode  II  it  can  be  seen  that  the  crack-tip  opening 
decreases  with  increasing  mode  II  contribution  and  the  crack  tip  tends  to  rotate  in  a 


clockwise  direction,  which  reaches  its  extremum  in  the  case  of  pure  mode  II.  In  the  case  of 
high  mode  II  values,  the  extreme  hydrostatic  state  of  stress  around  the  tip  cannot  be  relieved 
by  crack-tip  blunting  (as  in  cases  with  overwhelming  mode  I  contributions)  where  finite 
strains  create  a  relatively  sm^th  crack  tip. 

Figures  49  to  52  depict  only  the  upper  and  lower  element  layers  around  the  crack  tip  in 
deformed  and  undeformed  states  (with  a  magnification  factor  of  10)  for  selected  cases  of 
mixed-mode  loading.  The  influence  of  the  rotation  of  the  crack  tip  and  the  deviation  from  the 
center  line  of  the  undeformed  crack  is  distinct  for  cases  of  high  mode  II  values. 

Figures  53  to  56  show  deformed  versus  undeformed  meshes  of  the  outer  region  between 
1  and  10  mm  radii  around  the  crack. 


The  strain  energy  density  criterion,  according  to  Sih  [31,32],  has  not  only  been 
capable  of  predicting  fracture  under  brittle  material  behavior,  but  also  fracture  in  the 
elastic-plastic  regime. 

The  strain  energy  density,  using  the  notation  of  Sih,  is  given  as  : 


/°W\  f  ' 
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(9.3) 


where  €jj  is  the  strain  tensor  and 
crjj  is  the  stress  tensor. 

The  fact  that  excessive  change  in  shape  can  be  associated  with  yielding  while  excessive 
change  in  volume  can  be  associated  with  fracture  lead  to  the  formulation  of  the  strain  energy 
density  criterion  for  elastic  plastic  material  behavior.  It  is  postulated,  [31,32],  that 


46 


-'""N  av; 


maximum  yielding  occurs  when  the  strain  energy  density  reaches  its  maximum  value 
(dW/dV)max  whereas  fracture  initiation  is  associated  with  minimum  strain  energy  density 

(dW/dV)min.  Chow  and  Xu  [31]  investigated  the  extension  of  the  modification  of  the  strain 
energy  criterion  in  the  elastic-plastic  regime  for  the  case  of  mixed-mode  loading.  It  was 
observed  that  the  strain  energy  density  criterion  led  to  incorrect  predictions  of  the  angle  of 
fracture  initiation  since  two  local  minima  can  be  observed  for  some  mixed-mode  cases. 

Figure  57  shows  the  variation  of  the  strain  energy  density  versus  the  angle  9  along  a 
circle  of  radius  r=0.4  mm  around  the  crack  for  selected  cases  of  mixed-mode  loading.  The 
minimum  strain  energy  density  is  located  at  angles  from  180°  for  pure  mode  I  to  90°  for 
pure  mode  II.  This  contradicts  experimental  evidence  [68-70].  Therefore  the  assumption  of 
the  crack  growth  direction  was  restated  by  Sih  [71]  that  the  direction  of  maximum  value  of 

(dW/dV)mln  governs  the  onset  of  crack  growth.  These  values  range  between  9  =0°  for  mode 

I  to  0  =»  -90°  for  the  case  of  of  pure  mode  II  and  are  given  in  Table  6  for  all  cases  of 
mixed-mode  loading  considered. 

In  contrast  to  the  maximum  values  of  the  strain  energy  density,  which  increase  sharply 
for  growing  mode  II  contribution  due  to  higher  stress  and  strain  components  around  the 
crack,  the  maximum  values  of  ( dW/dV)min  remain  for  all  cases  remarkably  constant. 

Very  good  agreement  between  the  location  of  (dW/dV)max  and  the  maximum  yield  can  be 
found  by  comparing  the  the  angular  position  of  the  highest  effective  stress  in  Figure  17  and 
Figures  20  to  26  with  the  predicted  values  by  the  strain  energy  density  criterion. 

The  singular  behavior  of  the  strain  energy  density  versus  the  distances  ahead  of  the 
cracK  along  the  line  9-0°  is  depicted  in  Figure  58  in  full  logarithmic  representation.  In 
equation  (8.5)  it  has  been  assumed  that  the  strain  energy  density  is  governed  by  a  1/r 
singualrity  inside  the  HRR  field.  It  is  evident  that  low  contributions  of  mode  II  values  result 
in  weaker  singularities  in  the  strain  energy  density  as  the  intermediate  zone  is  approached. 
Determination  of  the  powers  of  the  singularities  which  are  given  in  Table  7  demonstrate  that 
the  determined  powers  are  in  excellent  agreement  with  the  predicted  value  even  for 
intermediate  cases  of  mixed-mode  loading.  Cases  distinguished  by  low  mode  II  contribution 
show  clearly  the  weak  singular  behavior  inside  the  intermediate  zone.  As  the  crack  tip  is 
approacned,  it  can  be  seen  that  the  singular  behavior  of  these  cases  converges  towards  1/r. 
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This  observation  reinforces  the  correctness  of  the  assumption  in  equation  (8.5)  that  the 
HRR  field  solution  is  valid  if  the  strain  energy  density  exposes  a  1/r  singularity  inside  the 
plastic  zone  around  a  crack. 


The  stress  functions  can  be  resolved  from  near  tip  field  stress  and  displacement 
components  obtained  from  the  finite  element  results.  All  parameter  required  for  65  equally 
spaced  nodes  along  a  circular  path  of  radius  r  =  0.4  mm  around  a  crack  were  read  from  the 
ABAQUS  data  output  file.  Figure  18  shows  that  the  selected  radius  is  still  within  the 
dominant  crack-tip  singular  field  for  all  mixed-modes  cases.  Paths  closer  to  the  crack  tip 
showed  significant  scatter  in  the  values  of  <t„  and  c0g  components. 

Determination  of  the  stress  function  results  in  solving  the  three  equations  given  in  equation 
(8.8),  where  the  components  of  cr jj  and  Uj  are  given  in  polar  coordinates  : 


^ij  - - -  )  n+1  cr. . 
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The  equivalent  stress  cre  for  the  case  of  negligible  elasticity  in  the  case  of  plane  strain  is 


given  as  : 


5-e  =  V  3  (  o-rr-  5r00)  <  0-n 


(10.2) 


Due  to  the  dependence  of  ln  on  values  of  a,;  ,  Uj  and  its  angular  derivatives  (along  a 
circular  path  taken  from  -n  to  n),  the  solution  of  this  system  of  equations  can  only  be 
accomplished  iteratively.  In  the  present  case,  use  of  the  interval  halving  method  was  made 
which  showed  rapid  convergence  for  all  cases  considered.  Figure  59  explains  the  principle 
of  the  interval  halving  method. 

For  an  initial  estimate  of  ln  the  stress  and  displacement  functions  are  calculated 
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according  to  equation  (10.1).  The  derivatives  of  the  nodal  displacements  with  respect  to  0 
are  obtained  using  a  sixth  order  finite  difference  formula  and  the  integration  is  performed 
using  Simpson's  second  order  integration  rule.  For  incrementally  increasing  values  of  ln. 
both  stress  and  displacement  components  and  derivatives  are  calculated  and  substituted  into 
equation  (10.2)  to  gain  a  new  value  of  ln.  The  calculated  value  of  ln  approaches  the  assumed 
value  of  ln  and  finally  surpasses  it.  Since  this  can  be  expressed  as  the  difference  (F) 
between  the  assumed  value  lna  and  the  calculated  value  lnc 

F  =  ^na  ‘  ^nc  •  (10-3) 

a  change  of  sign  in  (F)  is  expected  between  two  incremental  values  of  lna.  The  following 
interval  halving  procedure  narrows  the  interval  where  the  change  in  sign  has  occurred  down 
to  a  specified  residual.  A  numerically  accurate  value  for  ln  can  therefore  be  obtained. 

Figures  60  to  69  show  the  stress  functions  <rrr,cree  ,  crr0  and  the  effective  stress 
function  crQ  which  are  normalized  by  setting  the  maximum  value  of  the  0  -variation  of  the 
effective  stress  to  unity. 

For  all  cases  considered,  the  components  cr  ee  and  crr0  show  the  expected  value  of  zero 
on  the  crack  surfaces,  that  is,  for  angles  of  0  =■  ±  n.  In  contrast  to  pure  mode  I  in  the  linear 
elastic  case,  <xrr  assumes  positive  values  on  either  side  of  the  crack  flank.  On  opposite 
surfaces  of  the  crack  crrr  is  positive  for  both  pure  mode  I  and  a  mixed-mode  case 
distinguished  by  the  stress  intensity  ratio  of  K|/K u  ■  2222/178.  This  observation 
contradicts  Shih’s  statement  [27]  that  "  for  any  deviation  of  mode  I,  the  minus  sign  holds  in: 

cr rr  (0=n)*  -crrr  (0»-n)  (10.4) 

for  0  S  M*3  <  rc ." 

Also,  the  equivalence  of  the  magnitudes  of  radial  stress  components  on  the  surfaces  of  the 
crack  at  equal  distance  from  the  crack  could  not  be  shown  except  for  the  symmetic  cases  of 
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both  pure  mode  I  and  mode  II.  This  can  also  be  seen  from  the  magnitude  of  the  effective  von 
Mises  stresses  on  the  crack  surfaces  in  Figures  20  to  26.  Since  the  effective  stress  contains 
<jrr  as  the  only  component  different  from  zero,  their  contours  had  to  match  up  at  the  crack 

flanks.  Even  though  the  observation  of  equal  radial  stresses  on  crack  surfaces  by  Rice  and 
Sudiansky  [71]  seems  plausible,  no  reason  for  this  mismatch  could  be  found.  This 
phenomenon  could  also  be  observed  in  other  studies  on  mixed-mode  fracture  [27,72,73]. 

From  equation  (8.7)  the  mixity  parameter  MPI  can  be  resolved  from  the  stress 
components  cr00  and  crr0  of  the  finite  element  results.  Some  scatter  in  the  ratios  of  cr00 

over  crr0  for  distances  of  less  than  0.5  mm  from  the  crack  tip  was  observed  and,  therefore, 

a  least  square  approximation  was  employed  through  ail  datapoints  for  distances  ranging  from 
0.08  to  1  mm  from  the  crack  tip. 

Figure  70  shows  the  ratios  of  cr00  over  <xre  of  same  integration  points  along  the  crack 

and  their  corresponding  least  square  approximations  for  three  mixed-mode  cases.  The 
constant  character  of  the  curves  is  evident. 

The  relationship  between  the  calculated  values  of  the  mixity  parameter  MPI  and  the 
constant  ln  given  in  equation  (8.6)  is  Shown  in  Figure  71.  The  results  compare  well  with 
those  generated  by  Shih  [27]  for  the  hardening  powers  of  n  -  5  and  n  -  13. 
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1L  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 


In  the  present  study,  local  crack-tip  quantities  were  investigated  for  general 
mixed-mode  crack  problems  under  the  condition  of  plane  strain.  In  order  to  become 
independent  of  a  specific  specimen  geometry,  the  local  crack-tip  region  was  modeled  as  a 
disk  with  the  crack  tip  at  its  center.  Displacements  on  the  boundary  of  the  local  crack 
region  were  calculated  from  assumed  combinations  of  stress  intensity  factors  for  mode  I 
and  mode  II. 

The  strain  energy  density  criterion,  according  to  Sih,  was  applied  as  a  fracture 
criterion.  This  provides  a  concise  relationship  between  a  given  strain  energy  density 
factor  and  the  stress  intensity  factors,  K|  and  K||.  For  ten  comparable  cases  of  loading 
(which  span  the  range  from  pure  mode  !  to  pure  mode  II)  the  body  response  for 
elastic-plastic  material  behavior  was  calculated  using  the  finite  element  package  ABAQUS. 
The  eight  node  plane  strain  isoparametric  elements  employed  were  degenerated  into 
triangular  shaped  elements  around  the  crack  tip  to  produce  a  1/r  singularity  in  strains. 
The  material  data  of  the  stainless  steel  A304  was  used  in  the  finite  element  calculation. 

The  path  independent  J1 -integral  (which  is  a  governing  parameter  of  the  amplitude  of 


the  crack-tip  singularities  of  the  stress  and  strain  fields)  was  calculated  according  to  the 
method  of  virtual  crack-tip  extension  that  is  available  in  ABAQUS.  The  direct  integration 
method  was  applied  for  four  cases  of  mixed-mode  loading,  which  had  to  be  subdivided  into 
separate  consecutive  steps  to  reach  convergence.  For  these  cases  the  virtual  crack 
extension  method  could  not  be  applied  by  ABAQUS.  The  J1 -integral  value  was  calculated 


using  the  direct  integration  method  along  nine  circular  paths  whose  radii  from  the  crack 
tip  spanned  between  0.6  and  8.3  mm.  Generally,  very  good  path  independence  was  observed 
which  indicated  the  correctness  of  the  obtained  values  and  justified  for  their  further  use. 
Good  agreement  between  the  J1 -integral  values  obtained  by  both  methods  was  observed 


when  both  were  calculated.  Deviations  between  the  results  obtained  by  both  methods  were 
within  5.7  percent  for  all  cases  considered. 

J 2-integral  values,  which  have  attained  limited  usage  in  the  field  of  fracture 

mechanics,  were  also  evaluated.  Less  path  independence  especially  for  cases  of  more 
balanced  mode  I  and  mode  II  contributions  was  observed  due  to  larger  numerical  errors  in 
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the  evaluation  of  the  displacements  which  was  caused  by  higher  displacement  gradients  in 
the  y  -direction. 

The  characteristic  singular  behavior  of  stresses  and  strains  within  the  HRR  field,  the 
intermediate  zone,  and  the  far  field  was  best  expressed  by  the  dependence  of  the  effective 
von  Mises  stress  on  the  distance  from  the  crack  tip,  along  the  line  0-  0°.  Extraction  of  the 
powers  characterizing  the  singular  behavior  of  stresses  and  strains  indicated  the 
dependence  of  the  size  of  the  HRR  field  on  increasing  mode  II  contributions.  In  the  case  of 
pure  mode  I,  the  characteristic  exponent  of  -1/7  (for  the  given  material)  in  the  effective 
von  Mises  stress  was  never  reached  due  to  the  distinct  influence  of  crack-tip  blunting. 
Extreme  hydrostatic  stresses  in  the  vicinity  of  the  crack  tip  are  reduced  by  finite  strains 
and  the  stress-free  crack-tip  surfaces.  For  higher  mode  II  contributions,  the  crack  width 
decreased  significantly  to  an  almost  sharp  crack  in  the  case  of  pure  mode  II.  Larger  plastic 
zone  sizes  and  less  crack-tip  blunting  result,  therefore,  in  a  distinct  HRR-field  of 

increasing  size  which  was  observed  to  be  valid  for  a  distance  of  approximately  25  J/<ry  in 
the  case  of  pure  mode  II.  Accordingly,  the  effective  stress  increased  sharply  for  higher 
mode  I!  contributions. 

The  strain  energy  density  was  investigated  to  determine  its  applicability  as  a  fracture 
criterion  for  elastic-plastic  material  behavior  under  mixed-mode  loading.  It  was  stated  by 
Sih  [29]  that  the  maximum  yield  occurs  where  the  strain  energy  reaches  its  maximum 
which  could  be  verified  in  this  study.  The  modified  formulation  of  the  condition  that  crack 
extension  occurs  under  the  angle  where  the  minimum  of  the  strain  energy  density  shows  a 
maximum  seems  to  be  promising.  For  increasing  mode  II  values  the  direction  of  crack 
growth  initiation  could  be  shown  to  move  from  0°  to  -90°  relative  to  the  crack  plane. 
This  is  comparable  to  the  results  of  the  linear  elastic  solution. 

An  emphasis  of  the  present  study  was  the  numerical  extraction  of  the  stress  functions 
from  the  finite  element  solution,  for  a  given  material  with  a  hardening  factor  of  six.  The 
stress  functions  obtained  compared  reasonably  well  with  those  obtained  by  Shih. 
Theoretical  considerations  suggest  radial  stresses  of  equal  magnitude  for  equal  distances 
from  the  crack  tip  on  either  crack  surface.  This  would  imply  that  the  radial  stress 
function  (which  is  the  only  component  different  from  zero  on  the  crack  surface)  has  the 
same  magnitude  at  either  side  of  the  crack  surface.  This  could  only  be  shown  for  ihe 
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symmetric  mode  I  case  and  the  skew-symmetric  mode  II  case.  In  all  other  cases  of  mixed¬ 
mode  loading  a  deviation  in  the  magnitude  of  the  radial  stress  component  on  either  side  of 
the  crack  could  be  observed.  The  effective  von  Mises  stress  (which  is  only  dependent  on  the 
radial  stress  along  the  surfaces  of  the  crack)  showed  this  mismatch  also.  For  any  amount 
of  mode  II  contribution,  the  predicted  change  in  sign  of  the  radial  stress  acting  on  opposite 
crack  surfaces  was  not  verified  by  the  finite  element  results.  One  case  investigated,  where 
the  specimen  was  subjected  to  a  load  resulting  in  stress  concentration  factors  of  K|  *  2222 
and  K||  -  1 78  NVmm3,  showed  positive  radial  stresses  on  either  side  of  the  crack  surface. 
The  finite  element  calculations  indicated  that  the  sign  of  crrr  on  the  crack  surface  jumps 

suddenly  from  positive  to  negative  values  for  a  more  intermediate  mixed-mode 
combination. 

The  investigation  of  the  observed  discrepancy  between  the  behavior  of  theoretical 
radial  stresses  along  the  crack  surface  and  the  finite  element  solution  was  not  studied  in 
this  thesis.  An  important  extension  of  this  work  should  include  research  on  this 
phenomenon. 

An  extension  of  this  study  should  include  investigations  of  the  local  crack-tip 
quantities  of  a  suggested  mixed-mode  specimen.  The  mixed-mode  fracture  specimen  (due 
to  Richard  [34])  has  been  shown  to  simulate  very  well  arbitrary  mixed  modes  ranging 
from  pure  mode  I  to  pure  mode  II  and  should  be  given  consideration.  Both  numerical  and 
experimental  studies  are  necessary  to  further  investigate  two  interesting  concepts  of 
fracture  criteria  which  predict  both  the  onset  of  crack  growth  and  the  direction  of  crack 
extension  for  ductile  materials  under  mixed-mode  loading  : 

i )  The  strain  energy  density  criterion  according  to  Sih 

The  modified  formulation  of  the  predicted  angle  of  fracture  initiation  appears  to  be 
promising  in  the  finite  element  results  but  needs  further  experimental  confirmation. 
Further  investigation  of  crack  growth  initiation  in  relation  to  its  associated 
strain  energy  density  ,  especially  in  the  range  of  high  mode  II  values  seems  nece  ;sary. 

i  i )  The  T-criterion  according  to  Theocaris 

This  criterion  (which  has  been  extended  very  recently  for  mixed-mode  loacings 


54 


k » .  • 


f, v* .  * 


51 

js 

fi  j 

s! 

r , 

i«d 


under  ductile  material  behavior)  incorporates  the  HRR  field  solution.  Experiments 
with  inclined  cracks  in  a  center  cracked  panel  showed  good  agreement  in  the  prediction 
of  the  angle  of  extension  versus  the  crack  and  failure  loads.  Both  finite  element 
analysis  and  experimental  investigations  with  a  specimen  which  can  reproduce  all 
combinations  of  mode  I  and  II  may  establish  this  criterion  in  the  field  of  fracture 
mechanics. 
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FiQura  3  : 


Definition  of  crack  angle  and  fracture  angle  in  thj  center  cracked  panel 
with  slanted  crack  under  umaxi&i  tensile  stress. 
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Figure  6  :  Mapping  of  the  eight-node  prabolic  element  from  spacial  coordinates 

(x,  y)  to  local  coordiantes  (5.  n). 


“l*  •  »  *  »  „  »  „  »  ,  »  -  »  W  •'.k  ’  »  >  * 

*  **— *■  AA  *4.-^  ‘  ‘  A+i  j4-.\  JL./^  M- A  A /i  ^  ^ 


•  m  a  a-i  - *  •>: -j^u  A 1 


Fiflurfl  B  •  Generation  of  the  eight-node  collapsed  element  : 

( a )  rectangular  eight-node  element, 

( b  )  degenerated  triangular  eight-node  element. 
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0/2018 


Degrees  (X  10**2) 

Effective  von  Mises  stresses  around  crack  tip.  Radius  r  «  0.4  mm  for 
selected  mixed  mode  cases. 


for  all  of  mixed  mode  cases  considered. 
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K|/K„  -  2230/0 


ratio  :  K,/K„  -  1927/987 


Equivalent  von  Mises  stress  around  crack  tip.  fl 
Stress  intensity  factor  ratio  :  K/K|(  -  1405/1462 


Equivalent  von  Mises  stress  around  crack  tip. 
Stress  intensity  factor  ratio  :  K/Kf|  =  772/1774 


crack  tip.  Radius  r=  1  mm. 


Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm.  outer 


Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 


Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm.  outer 


Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm. 


Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm. 

Stress  intensity  factor  ratio  :  K/K(|  =  396/1903. 


Equivalent  von  Mises  stress  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm. 


Equivalent  plastic  strain  around  crack  tip.  Radius  r«  1  mm. 
Stress  intensity  factor  ratio  :  K/K(|  -  1927/987. 


Equivalent  plastic  strain  around  crack  tip.  Radiu: 
Stress  intensity  factor  ratio  :  K/K||  -  1098/1633. 


around  crack 


Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K,/K||  =  1927/987. 
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Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm.  outer 
radius  10  mm.  Stress  intensity  factor  ratio  :  K/Klt  =  772/1774. 


Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm.  Stress  intensity  (actor  ratio  :  K(/Kj|  -  396/1903. 
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Equivalent  plastic  strain  around  crack  tip.  Inner  radius  1  mm,  outer 
radius  20  mm.  Stress  intensity  factor  ratio  :  K|/K||  *=  0/2018. 
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Deformed  Mesh,  outer  radius  r»  1mm. 
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Deformed  Mesh,  outer  radius  r«  1mm. 


Deformed 


deformation 


Crack  deformation.  Plotted  crack  length  10  mm. 
Stress  intensity  factor  ratio  :  K/K||  -  1098/1633. 


|/K|.  =  1098/1633. 


variation  ot  the  stress  lunctions  a 


variation  of  the  stress  functions  cr 


selected  mixed-mode  cases. 


No. 

K|  -  Value 
[N/Vmm3] 

K||  -  Value 
[N/Vmm3] 

Mixity  Parameter 
(Me) 

1 

2230.0 

0. 

1. 

2 

2222.9 

175.8 

0.949 

3 

2107.7 

670.8 

0.804 

4 

1927.0 

987.6 

0.698 

5 

1683.8 

1252.0 

0.584 

6 

1405.6 

1462.2 

0.487 

7 

1098.7 

1633.8 

0.377 

8 

772.2 

1774.5 

0.261 

9 

396.3 

1903.8 

0.130 

10 

0. 

2018.3 

0. 

pure  mode  I 


pure  mode  II 


Table  1  : 


K|  -  K||  -values  according  to  the  fracture  criterion  by  Sih.  The 
elastic  mixity  parameter  Me  is  given  by  equation  (8.1) 


Ni  *  7-  0  -  S  )  (1  -  *n  )  (  -1  -  H  -  n  ) 
4 


n2  *  7-  O  +  5  )  (1  -  r>  )  (  -1  -+?  -  n  ) 
4 


n3  *  —  0  +  S  )  (1  +  n  )  (  -1  +  S  +  ri  ) 


N4  ■  —  (1  -  ?  )  (1  +  n  )  {  -1  -  S  +  r»  ) 


NS  -  1  (i  -  ?2)  (i  -  ^  ) 


n6  ,  1  (i  +  ?  )  (i  -  ^2  ) 


N7-  j(1-52)(1  +ti  ) 


n8  -  1  (i  -  e  )  (i  -  n2 ) 


Table  2  : 


Shape  functions  of  the  eight-node  isoparametric  finite  element. 
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Chemical  Composition  in  Weight  Percent : 

Element:  C  Mn  P  S  Cu  Si  Ni  Cr  Mo  Co  N 
0.03  1.62  0.034  0.015  0.31  0.43  10.04  18.3  0.31  0.15  0.07 


Material  Data  : 


Property  : 

Symbol 

Unit 

Value 

Young's  modulus 

:  E 

MPa 

203800. 

Yield  stress 

:  <ry 

MPa 

265. 

Maximal  stress 

:  ^max 

MPa 

640. 

Poisson's  ratio 

:  vj 

- 

0.3 

Critical  stress  intensity  factor 

:  KIC 

N/Vmm3 

4300. 

Critical  J-integral  value 

:  J1C 

N/mm 

108. 

TafrlQ  2  : 


Chemical  composition  and  material  data  of  the  stainless  steel  A304. 


J-Integral  Value 
(Vitual  Crack  Extension) 

[N/mm] 


J-integral  Value 
(Direct  Integration) 

[N/mm] 


Deviation 


% 


20.49 

20.50 

0.05 

20.97 

21.17 

0.95 

- 

20.21 

- 

- 

18.94 

- 

- 

17.19 

- 

- 

15.46 

- 

11.57 

11.88 

2.61 

9.94 

10.22 

2.74 

9.28 

9.87 

5.77 

9.19 

9.73 

5.55 

Comparison  oi  J-integral  values  obtained  by  the  virtual  crack 
extension  method  and  the  direct  integration  method.  ' 


Mixed  Modes  : 


Number : 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

K|: 

2230.  2223. 

2107.  1927. 

1638. 

1405. 

1098. 

772. 

396. 

0. 

K|| : 

0. 

176. 

671. 

988. 

1252. 

1462. 

1633. 

1775. 

1904.  2018. 

HRR-Zone  : 

inner  radius  : 

0.08 

0.08 

0.08 

0.08 

0.1 

0.1 

0.1 

0.1 

0.2 

0.2 

outer  radius : 

0.12 

0.12 

0.3 

1. 

4. 

6. 

8. 

8. 

9. 

9. 

power : 

-0.114 

-0.119 

-0.134 

-0.134 

-0.135 

-0.140  -0.140 

-0.141 

-0.142 

-0.143 

lolarm-Zang 

inner  radius  : 

0.15 

0.15 

0.9 

2. 

5. 

8. 

outer  radius  : 

1. 

1. 

2.2 

5. 

10. 

12. 

- 

- 

- 

- 

power : 

-0.0391 

-0.041 

-0.041 

-0.060 

-0.041 

-0.053  - 

- 

- 

- 

E'astic  Zone: 

inner  radius  : 

2. 

5. 

6. 

- 

- 

- 

- 

- 

- 

. 

outer  radius  : 

12. 

12. 

12. 

- 

- 

- 

- 

- 

- 

. 

power : 

-0.51 

-0.499 

-0.50 

. 

. 

Table  5  :  Powers  characterizing  singular  behavior  of  the  effective  von  Mises 

stress  along  the  line  9  =  0°  ahead  of  the  crack  tip.  (  All  distances  are 
given  in  mm  ) 


y.-.'V 


Mixed  Mode  Cases  : 


No. 

K| 

Kll 

Predicted  Fracture 

Magnitude  of  Strain 

[N/Vmm3] 

Angle 

Energy  Density 
[MJ/m3] 

1 

2230. 

0. 

0. 

1.83 

2 

2222.9 

175.8 

-10.48 

1.82 

3 

2107.7 

670.8 

-21.60 

1.71 

4 

1927.0 

987.6 

-33.39 

2.00 

5 

1638.8 

1252.0 

-45.00 

2.22 

6 

1405.6 

1462.2 

-50.80 

2.47 

7 

1098.7 

1633.8 

-68.03 

1.89 

8 

772.2 

1774.5 

-79.15 

1.95 

9 

396.3 

1903.8 

-84.56 

1.82 

10 

0. 

2018.3 

-90.00 

1 .925 

Table  6  :  Fracture  angle  0  and  corresponding  strain  energy  density  0.4 
mm  from  the  crack  tip  for  all  mixed-mode  cases  considered. 


Mixed  Mode  Cases : 


No. 

K| 

K|| 

Range  of  Distance  from 

Power  of  the 

Correlatit 

[N/Vmm3] 

Crack  Tip  Considered  [mm] 

Singularity 

1 

2230.0 

0. 

0.04- 

0.12 

-0.622 

-0.9999 

2 

2222.9 

175.8 

0.04  - 

0.12 

-0.699 

-0.9996 

3 

2107.7 

670.8 

0.04- 

0.20 

-0.880 

-0.9994 

4 

1927.0 

987.6 

0.04- 

0.40 

-0.943 

-0.9997 

5 

1638.8 

1252.0 

0.04- 

0.80 

-0.983 

-0.9998 

6 

1405.6 

1462.2 

0.04- 

1.0 

-1.030 

-0.9999 

7 

1098.7 

1633.8 

0.04- 

1.0 

-1.029 

-0.9999 

8 

772.2 

1774.5 

0.04- 

1.0 

-1.030 

-0.9998 

9 

396.3 

1903.8 

0.04- 

1.0 

-1.002 

-0.9998 

10 

0. 

2018.3 

0.04- 

1.0 

-1.000 

-0.9999 

Iahlfi-Z  :  Powers  of  the  singularity  of  the  strain  energy  density  in  the  vicinity 

of  the  crack  tip  along  the  line  0=0°  from  least  square  approximation. 
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In  the  present  work,  slow  stable  crack  growth  in  a  ductile  material  (A533B  steel)  is 
simulated  numerically  with  a  widely  used  state  of  the  art  commercial  finite  element 
code  (ABAQUS).  The  finite  element  formulation  uses  the  J2  flow  rule  of  incremental 
plasticity  based  on  small  deformation  theory.  An  experimentally  obtained  load  versus 
crack  growth  relation  is  employed  as  input  and  the  finite  element  mesh  models  the 
upper  half  of  a  modified  compact  tension  specimen.  The  material  behavior  is  modeled 
with  power  law  hardening  (n  -  10).  Two  meshes  are  used  for  the  crack  growth 
simulation:  a  coarse  mesh  with  1874  elements  and  1946  nodes,  and  a  fine  mesh  with 
2914  elements  and  3002  nodes.  The  finite  element  meshes  consist  of  four  -  node 
bilinear  plane  strain  isoparametric  elements  with  eight  degrees  of  freedom  (type 
CPE4).  Crack  extension  of  eight  millimeters  for  the  coarse  mesh  and  one  millimeter 
for  the  fine  mesh  were  simulated.  The  nodal  release  technique  is  used  as  a  numerical 
crack  growth  simulation  technique. 

The  evaluation  of  the  results  emphasizes  two  points:  (i)  to  determine  how  well  two 
crack  growth  criteria  (J/CTOD  and  CTOA)  characterize  mode  1  crack  extension  under 
plane  strain  conditions  and  in  what  interval  range  they  are  applicable,  and  (ii)  to 
examine  the  singular  field  variables  at  the  onset  of  crack  growth  and  associated  with 
the  quasi  -  static  crack  extension.  The  purpose  of  the  second  evaluation  is  to  determine 
whether  the  singularity  fields  suggested  by  asymptotic  methods  exist  independently 
from  each  other  (whether  a  transition  point  between  two  different  singularity  fields 
can  be  identified)  or  whether  a  superposition  of  the  singularity  fields  occurs  ahead  the 
stable  advancing  crack  tip. 


L. 


7.1  J  -  INTEGRAL  INTRODUCTION 

7.1.1  J  -  INTEGRAL  DETERMINED  BY  THE  VIRTUAL  CRACK  EXTENSION 
METHOO 

7.1.2  J  -  INTEGRAL  AS  CRITERION  FOR  STABLE  CRACK  GROWTH 

7.1 .3  CALCULATION  OF  THE  J  -  INTEGRAL 

7.1 .4  DISCUSSION  OF  THE  J  -  INTEGRAL  AS  CRITERION  FOR  STABLE  CRACK 
GROWTH 

7.2  CRACK  PROFILE  GEOMETRY 

7.2.1  DEFINITION  OF  THE  CRACK  PROFILE  PARAMETER 

7.2.2  DISCUSSION  OF  THE  CRACK  PROFILE  PARAMETER 

7.2.3  CRACK  TIP  OPENING  DISPLACEMENT  (CTOD)  -  RESULTS 
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Figure  1 :  The  three  local  independent  kinematic  movements  of  the  upper  and  lower 
surface  with  respect  to  each  other. 

Figure  2:  Idealized  constitutive  material  behavior. 

Figure  3:  Cottrell's  model  of  cleavage  fracture. 

Figure  4:  The  schematic  representation  of  a  commonly  used  simple  microscopic 
fracture  criterion. 

Figure  5:  Local  crack  growth  simulation  in  three  dimensions. 

Figure  6:  Stiffness  reduction  accomplished  with  spring  elements. 

Figure  7:  Load  versus  crack  growth. 

Figure  8:  Step  -  wise  crack  growth  simulation. 

Figure  9:  Stress  -  strain  curve  {A533B  steel). 

Figure  1 0:  Dimensions  of  the  employed  Compact  Tension  specimen. 

Figure  11a  -  lib:  Finite  element  modeling  of  the  Compact  Tension  specimen 
{coarse  mesh). 

Figure  12a  -  12c:  Finite  element  modeling  of  the  Compact  Tension  specimen 
(fine  mesh). 

Figure  13:  Boundary  conditions. 

Figure  14:  Eight  d.  o.  f.  bilinear  element. 

Figure  15:  Mapping  of  the  eight  d.  o.  f.  bilinear  element  from  natural  coordinates  into 
xy  -  space. 

Figure  16:  Locations  of  the  Gauss  quadrature  integration  points  within  an  eight  d.  o.  I. 
bilinear  isoparametric  element. 

Figure  17:  Definition  of  the  line  of  the  J  -  integral  around  a  notch. 

Figure  18:  Typical  J  •  resistance  curve. 

Figure  19:  Schematic  of  the  fields  surrounding  a  growing  crack  [34], 


Figure  20:  Evaluation  line  for  the  J  -  integral. 

Figure  21:  Verification  paths  for  the  path  independence  of  the  J  -  integral. 

Figure  22:  J  -  integral  versus  crack  growth. 

Figure  23:  Load  -  line  displacement  versus  crack  growth. 

Figure  24:  Tangent”  -  definition  of  the  crack  tip  opening  displacement. 

Figure  25:  Crack  profile  (coarse  mesh,  0  -  1  mm  crack  growth). 

Figure  26:  Crack  profile  (coarse  mesh,  1.25-  3  mm  crack  growth). 

Figure  27:  Crack  profile  (coarse  mesh,  4  -  8  mm  crack  growth). 

Figure  28:  Crack  profile  (fine  mesh,  0  -  1  mm  crack  growth). 

Figure  29a  to  29q:  Crack  tip  region  shown  for  the  advancing  crack  (coarse  mesh, 
magnification  factor  1). 

Figure  30a  to  30i:  Crack  tip  region  shown  for  the  advancing  crack  (fine  mesh, 
magnification  factor  1). 

Figure  31 :  Crack  tip  opening  displacement  versus  crack  growth  (coarse  mesh). 

Figure  32:  Crack  tip  opening  displacement  versus  crack  growth  (fine  mesh). 

Figure  33:  Crack  tip  opening  angle  versus  crack  growth  (coarse  mesh). 

Figure  34:  Schematic  representation  of  the  sensitivity  of  the  CTOA  dependent  on  the 
element  size  at  the  crack  tip. 

Figure  35:  Crack  tip  opening  angle  versus  crack  growth  (fine  mesh). 

Figure  36:  "45°*  definition  of  the  crack  opening  displacement. 

Figure  37:  Dependence  of  dn  on  n  and  cr<yE  for  plane  strain  [18]. 

Figure  38:  Prandtl  slipline  fields  for  a  steadily  growing  crack. 

Figure  39:  Cartesian  stress  components  at  a  crack  tip. 

Figure  40:  Strain  in  y  -  direction  along  the  line  0-0°  for  the  advancing  crack  (fine 
mesh,  0  -  1  mm  crack  extension). 
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Figure  41: 

Figure  42: 

Figure  43: 

Figure  44: 

Figure  45: 

Figure  46: 

Figure  47: 

Figure  48: 

Figure  49: 

Figure  50: 

Figure  51: 

Figure  52: 
Figure  53: 

Figure  54: 

Figure  55: 


Strain  in  y  -  direction  along  the  line  ©  »  0°  for  the  advancing  crack 
(coarse  mesh,  0  -  1  mm  crack  extension). 

Strain  in  y  -  direction  along  the  line  ©  ■  0°  for  the  advancing  crack 
(coarse  mesh,  1.5-4  mm  crack  extension). 

Strain  in  y  -  direction  along  the  line  ©  -  0°  for  the  advancing  crack 
(coarse  mesh,  5  -  8  mm  crack  extension). 

Stress  in  y  -  direction  along  the  line  ©  -  0°  for  the  advancing  crack 
(fine  mesh,  0  -  1  mm  crack  extension). 

Stress  in  y  -  direction  along  the  line  ©  -  0°  for  the  advancing  crack 
(coarse  mesh,  0  -  1  mm  crack  extension). 

Stress  in  y  -  direction  along  the  line  ©  -  0®  for  the  advancing  crack 
(coarse  mesh,  1.5-4  mm  crack  extension). 

Stress  in  y  -  direction  along  the  line  ©  -  0°  for  the  advancing  crack 
(coarse  mesh,  5  -  8  mm  crack  extension). 

Von  Mises  equivalent  stress  along  irne  ©  -  0°  for  the  advancing  crack 
(coarse  mesh,  0  -  1  mm  crack  extension). 

Von  Mises  equivalent  stress  along  line  ©  -  0°  for  the  advancing  crack 
(coarse  mesh,  1.5-4  mm  crack  extension). 

Von  Mises  equivalent  stress  along  line  ©  »  0°  for  the  advancing  crack 
(coarse  mesh,  5  -  8  mm  crack  extension). 

Von  Mises  equivalent  stress  along  line  ©  -  0°  for  the  advancing  crack 
(fine  mesh,  0  -  1  mm  crack  extension). 

Transition  point  position  versus  crack  growth  . 

Log  /  log  representation  of  the  von  Mises  equivalent  stress  along  line 
©  -  0°  for  the  advancing  crack  (coarse  mesh,  0  -  1  mm  crack 
extension). 

Log  /  log  representation  of  the  von  Mises;  equivalent  stress  along  line 
©  -  0°  for  the  advancing  crack  (coarse  mesh,  4  and  8  mm  crack 
extension). 

Elastic  strain  energy  density  along  line  ©  -  0°  -  180°  for  the  advancing 
crack  (fine  mesh,  0  -  .75  mm  crack  extension). 
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Figure  56:  Elastic  strain  energy  density  along  line  ©  -  0°  -  180°  for  the  advancing 
crack  (coarse  mesh,  1  •  3  mm  crack  extension). 

Figure  57:  Elastic  strain  energy  density  along  line  ©  -  0°  -  180°  for  the  advancing 
crack  (coarse  mesh,  4  -  8  mm  crack  extension). 

Figure  58:  Plastic  strain  energy  density  along  line  ©  -  08  -  180°  for  the  advancing 
crack  (coarse  mesh,  0  -  1  mm  crack  extension). 

Figure  59:  Plastic  strain  energy  density  along  line  ©  -  0°  -  180°  for  the  advancing 
crack  (coarse  mesh,  1.5-4  mm  crack  extension). 

Figure  60:  Plastic  strain  energy  density  along  line  ©  »  0°  -  180®  for  the  advancing 
crack  (coarse  mesh,  5  •  8  mm  crack  extension). 

Figure  61a  -  61m:lso  contours  of  the  von  Mises  equivalent  stress  for  the  advancing 
crack  (coarse  mesh). 

Figure  62:  Behavior  of  the  characteristic  radius  Rc  of  the  strain  intense  region 
(sharply  bordered  region  where  the  rate  of  energy  dissipation  is  high). 


Table  1 :  Relation  between  the  external  (applied)  load  and  the  crack  growth  (coarse 
mesh). 

T able  2:  Material  composition  of  A533B  steel  [66]  (in  weight  percent). 

Table  3:  Stress  -  strain  properties  of  A533B  steel. 

Table  4:  J  -  values  for  different  paths. 

Table  5:  J  -  integral  values  over  crack  extension. 

Table  6:  Comparison  of  the  results  for  the  J  -  integral  performed  in  this  work  with 
the  work  of  Hoff  [37]. 
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1.  INTRODUCTION 

Under  normal  circumstances,  a  structural  analysis  assumes  that  the  materials 
involved  are  ideal  homogeneous  flawless  materials,  (i.e.  stresses  and  strains  are 
uniformly  distributed  throughout  a  body).  Inglis  [1]  first  emphasized  the  significance 
of  intense  and  localized  concentrations  of  stresses  around  sharp  notches.  Neuber  [21 
resolved  this  observation  of  stress  concentrations  caused  by  notches  by  introducing 
stress  concentration  factors. 

Griffith  contributed  pioneering  work  to  this  subject  in  the  early  1920's  [3.41.  He 
developed  a  continuum  mechanics  based  formulation  of  the  change  in  strain  energy  due 
to  the  presence  of  a  crack  in  brittle  elastic  solid.  Often  this  work  is  quoted  as  the 
starting  point  of  fracture  mechanics  as  an  independent  branch  of  mechanics.  Sneddon 
[51  deduced  expressions  for  the  stress  distribution  in  the  neighborhood  of  a  crack  in  an 
elastic  solid  from  complex  stress  functions  developed  by  Westergaard  [6J. 

The  next  step  in  the  development  of  the  theory  of  fracture  mechanics  was  made  by 
Irwin  in  the  1950's.  He  observed  that  there  are  three  independent  local  kinematic 
movements  of  the  upper  and  lower  crack  surfaces  with  respect  to  each  other  (fig.l ) 
[71. 

1 )  Opening  Mode  or  Mode  1 

2)  Sliding  Mode  or  Mode  2 

3)  Tearing  Mode  or  Mode  3 

Essentially  all  stress  systems  in  the  near  crack  -  tip  region  may  be  derived  from  these 
three  modes  of  loading.  Since  the  opening  mode  (or  mode  1)  represents  the 
predominant  stress  situation  in  many  practical  cases,  most  of  the  research  is  done  in 
this  area.  Building  on  the  associated  stress  fields  in  the  near  crack  -  tip  region  of  the 
three  different  crack  movements,  Irwin  deduced  the  stress  intensity  factor  (K) 
concept  [8J,  where  K  describes  the  intensity  of  the  elastic  crack  -  tip  stress  field. 
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Previously  Orowan  recognized  that  for  relatively  ductile  materials,  the  work  done  in 
plastic  deformation  is  much  larger  then  the  energy  required  to  form  new  crack  surface 
[9].  From  these  observations  Irwin  defined  a  material  property  G  which  is  the  total 
energy  released  during  crack  extension  [10].  In  addition  he  demonstrated  the 
equivalence  of  G  and  K  for  linear  elastic  material  behavior.  This  property  is  the  basis 
of  brittle  fracture  mechanics  today.  Since  the  stress  distribution  characteristics 
around  a  crack  are  always  the  same,  material  properties  can  be  found  by  testing 
suitable  specimens.  Such  material  properties  like  Gc  (critical  energy)  or  Kc  (critical 
stress  intensity  factor),  once  found,  can  be  compared  with  the  G  or  K  value  of  a  body 
subjected  to  a  certain  load  condition  and  the  designed  structure  simply  has  to  satisfy 
the  following  conditions 


eq.  1.01 


K<^ 


eq.  1.02 


where  s  is  the  safety  factor.  Up  to  this  point,  linear  elastic  material  behavior  has  been 
used  as  basic  assumption  to  develop  the  theory.  Therefore  the  discipline,  using  this 
principle,  is  called  Linear  Elastic  Fracture  Mechanics  (LEFM). 


Wells,  in  the  early  1960's  [11],  introduced  the  concept  of  the  crack  opening 
displacement.  This  was  the  first  example  of  a  fracture  concept  developed  beyond 
general  yielding.  Wells's  work  provided  the  basis  for  the  semi  empirical  'COD  Design 
Curve'  approach,  used  today  (especially  in  the  United  Kingdom  )  for  fracture  under 
contained  yielding  conditions.  Hutchinson  [12]  and  Rice  and  Rosengren  [13]  derived 
(under  the  assumption  of  a  power  law  hardening  material  in  the  nonlinear  region  of 
material  behavior)  solutions  for  the  stresses  and  strains  near  a  crack  -  tip  using  the 
deformation  theory  of  plasticity.  Rice  subsequently  [14]  deduced  an  alternative  (but 


equivalent)  approach  to  the  COD  Design  Curve.  He  introduced  the  J  -  integral,  a  path 
independent  contour  integral  around  the  crack  -  tip.  Although  several  path  independent 
contour  integrals  have  been  advanced  independently  [15,  16,  17],  for  any  fracture 
mechanics  analysis  where  significant  plasticity  occurs,  either  the  COD  Design  Cun/e 
or  the  J  -  integral  is  used.  The  approach  essentially  is  the  same  as  in  linear  elastic 
materials  shown  in  eq.  1.01  and  eq.  1.02,  namely  the  designed  structure  has  to  satisfy 
the  condition 


The  discipline  of  fracture  mechanics  where  elastic  •  plastic  deformation  must  be  taken 
into  account  is  called  Elastic  Plastic  Fracture  Mechanics  (EPFM).  Due  to  the 
complexity  of  the  problems  in  EPFM,  progress  in  this  discipline  is  not  as  advanced  as 
in  LEFM. 

The  Finite  Element  Method  (FEM)  is  a  numerical  method  which  provides  the 
opportunity  to  simulate  elastic  and  plastic  material  behavior.  By  formulating  different 
types  of  idealized  constitutive  behavior  (not  only  nonlinear  elastic  corresponding  to 
deformation  theory  of  plasticity,  but  also  incremental  plastic  corresponding  to  a  flow 
theory  of  plasticity)  (fig.  2),  it  is  possible  to  characterize  a  fracture  within  a  body 
under  arbitrary  load  conditions.  With  this  tool,  bodies  subjected  to  complex  loading 
conditions  (in  elasticity  as  well  as  in  plasticity)  can  be  examined. 

In  fracture  mechanics  today,  engineering  calculations  are  not  limited  to  the 
determination  of  the  combination  of  the  critical  crack  size  -  load  conditions  for 
fracture  instability.  In  addition,  calculations  to  determine  the  rate  of  progression  of  a 
crack  are  performed.  There  are  several  distinct  types  of  crack  growth: 

-  fatigue  crack  growth, 

-  creep  crack  growth. 
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-  environmentally  assisted  crack  growth,  and 
•  stable  crack  growth. 


Fatigue  crack  growth  occurs  in  structures  which  are  operating  under  alternating  loads 
sufficiently  severe  to  make  fatigue  resistance  a  primary  design  criterion.  The 
approach  for  solving  this  problem  is  to  relate  the  change  in  crack  length  with  the 
number  of  applied  load  cycles.  A  widely  used  equation  for  this  relation  is  the  'Paris 
Law'  [18] 

-jJ-  -  g(AK)m  e<*-  104 

where  Kmaxis  the  upper  load  stress  intensity  factor, 

Kmjn  is  the  lower  load  stress  intensity  factor, 

AK  is  Kmax  -  Kmjn 

da  is  the  crack  length, 

dN  are  the  load  cycles,  and 

g,m  are  material  dependent  constants. 

One  difficulty  encountered  is  that  an  exact  definition  of  the  transition  from  initiation  to 
propagation  often  is  impossible. 

Creep  crack  growth  is  a  very  important  problem,  particularly  in  the  power  generating 
industry  and  aircraft  gas  turbines.  Metals  show  a  creep  behavior  at  temperatures 
greater  than  about  thirty  percent  of  their  absolute  melting  temperatures.  There  are 
two  competing  mechanisms  to  describe  the  time  dependent  crack  growth  behavior.  The 
first  mechanism  builds  upon  the  blunting  of  the  crack  -  tip.  This  phenomenon  is 
observed  experimentally  and  has  been  simulated  numerically.  Due  to  the  crack  -  tip 
blunting,  the  stress  field  ahead  of  the  crack  relaxes  and  tends  to  retard  crack  growth. 
The  other  mechanism  results  in  an  accumulation  of  creep  damage  in  the  form  of 
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microcracks  ahead  the  crack  -  tip.  These  microcracks  join  each  other  causing  the 
crack  to  extend. 

Environmentally  assisted  crack  growth  is  an  extremely  complex  problem  and  even 
experts  in  this  field  cannot  always  agree  on  the  precise  distinction  between  the 
different  types  of  environmental  cracking,  characterized  by  corrosion,  stress 
corrosion,  and  corrosion  fatigue.  Environmental  effects  on  fatigue  crack  growth 
strongly  depend  on  specific  material  •  environment  combinations,  as  well  as  on  the 
frequency  of  the  stress  cycle,  the  wave  form  of  the  stress  cycle  and  the  temperature.  In 
the  case  of  high  crack  growth  rates  the  environmental  effects  are  often  negligible. 

In  ductile  materials  (like  A533B  steel),  slow  stable  crack  growth  is  observed  after 
the  onset  of  crack  growth  due  to  extensive  plastic  deformation,  although  the  structure 
may  still  remain  in  service.  A  fracture  analysis  based  on  the  onset  of  crack  growth, 
therefore,  would  lead  to  an  overly  conservative  estimation  and  the  structure  would  be 
prematurely  removed  from  service.  The  problem  to  solve,  is  to  determine  what  amount 
of  stable  crack  growth  is  allowable  prior  to  the  onset  of  rapid  crack  propagation. 
Theoretical  foundations  for  this  subject  are  based  on  Elastic  Plastic  Fracture 
Mechanics  (EPFM).  Different  approaches  exist  to  solve  this  problem,  but  the  major 
obstacle  still  is  to  find  a  fracture  criterion  which  characterizes  stable  crack  growth 
after  crack  initiation. 

Parallel  to  the  macro  -  description  of  fracture  in  structures,  a  micro  -  mechanism 
approach  has  been  developed.  The  main  disadvantage  of  this  approach  is  the  lack  of 
experimental  verification  of  proposed  solutions.  Nevertheless,  much  research  has  been 
performed  in  the  past.  In  particular,  the  ability  to  relate  micro  -  mechanisms  of 
cleavage  and  ductile  fracture  to  the  fracture  mechanics  parameters  such  as  K,  J  - 
integral,  CTOD  and  CTOA  seems  to  be  the  key  for  a  successful  application  of  this 
concept. 


The  micro  -  mechanism  of  fracture  itself  is  divided  into  fast,  uncontrollable  crack 
extension  and  slow  stable  crack  growth.  Fast  crack  extension  occurs  below  the  cleavage 
transition  temperature.  This  cleavage  fracture  is  a  brittle  fracture  but  micro  - 
plasticity  is  not  excluded.  Transgranular  cleavage  fracture  occurs  in  structural  steels 
of  yield  strengths  generally  less  than  500  MPa,  while  intergranular  cleavage  occurs 
in  higher  strength  alloy  steels  [19].  Zener  [20]  suggested  that  there  is  an  array  of 
dislocations  at  the  initial  stage  of  crack  formation.  As  more  dislocations  enter  this 
array  they  are  squeezed  together,  producing  a  local  stress  concentration.  This  local 
stress  concentration  increases  until  a  crack  nucleus  is  generated.  Stroh  [21] 
presented  an  analysis  of  this  approach.  This  analysis  shows  that  cleavage  fracture 
would  not  be  predicted  using  this  dislocation  model.  Since  cleavage  fracture  is  observed 
experimentally,  the  model  proposed  by  Zener  appears  to  be  inadequate.  Cottrell  [23] 
suggested  a  mechanism  which  leads  to  easy  nucleation  in  bcc  metals.  In  this  rather 
straightforward  approach,  two  dislocations  intersect  on  the  cleavage  plane  and  form  a 
new  dislocation.  Equation  1 .05  describes  and  fig.  3  shows  this  mechanism. 

f  <TTl>(io,)tf<1,,>(lot,  -  a<001>(001,  ,05 

The  new  dislocation  has  a  lower  dislocation  energy  than  the  initial  one,  therefore, 
crack  nucleation  will  be  easy  and  crack  extension  is  explainable  by  connecting 
different  crack  nuclei. 

Above  the  fibrous/cleavage  transition  temperature  materials  behave  in  a  fully  ductile 
manner.  This  transition  temperature  for  A533B  steel  is  about  room  temperature. 
After  reaching  the  transition  temperature  the  crack  advances  by  the  coalescence  of 
voids.  These  voids  contain  inclusions  of  second  -  phase  as  well  as  nonmetallic  particles 
[24],  For  the  initiation  of  ductile  fracture,  a  simple  criterion  commonly  used  is  [25 
to  27] 

<5 jC  -  (0.5  to  2)  dp  eq.  1.06 


where  <j ic  is  the  critical  crack  opening  displacement  and 

dp  is  the  mean  void  initiation  particle  spacing. 

Figure  4  shows  this  microscopic  fracture  criterion  schematically.  Unlike  cleavage 
cracks,  this  ductile  material  behavior  is  based  on  cracks  which  are  too  blunt  to  be  able 
to  propagate  in  an  uncontrollable,  fast  way.  Local  microscopic  criteria  for  void  growth 
ahead  of  the  crack  -  tip  have  been  proposed  by  Green  and  Knott  [2d]  and  Rice  and 
Sorensen  [30]. 

In  the  present  work  a  macroscopic  crack  growth  study  is  performed.  Using  an  elastic  - 
plastic  (small  strain)  finite  element  analysis,  a  crack  in  a  compact  tension  specimen 
is  extended  quasi  statically  under  plane  strain  conditions.  The  material  employed  is  the 
bainitic  pressure  vessel  grade  steel  A533B  and  a  power  hardening  law  is  used  to 
represent  the  stress  •  strain  relationship.  The  von  Mises  equivalent  stress  is  used  as  a 
yield  criterion.  The  macroscopic  fracture  criteria  (J  /CTOD  and  CTOA)  are  examined 
as  to  their  usefulness  to  model  slow  stable  crack  growth. 

An  extended  evaluation  has  been  made  into  the  field  variables  in  the  vicinity  of  a  crack 
-  tip.  In  particular  the  changing  nature  of  the  field  variables  for  a  growing  crack  is 
examined  closely,  from  the  onset  of  crack  growth  to  eight  millimeters  of  crack 
extension. 


A  literature  survey  has  shown  that  three  different  finite  element  methods  are 
commonly  used  to  simulate  crack  growth: 

1)  Node  shifting, 

2)  Stiffness  reduction, 

3)  Nodal  release. 

Node  shifting  is  used  particularly  to  simulate  small  amounts  of  crack  growth,  namely 
less  than  one  element  size.  With  this  method  crack  blunting  can  be  simulated  very 
accurately  by  using  higher  order  elements.  If  larger  amounts  of  crack  growth  are 
needed,  it  is  possible  to  combine  node  shifting  with  nodal  release.  An  interesting 
application  is  the  simulation  of  local  three  dimensional  crack  growth  (fig.  5).  Neither 
the  nodal  release  nor  the  stiffness  reduction  methods  can  perform  this  simulation 
successfully. 

Stiffness  reduction  conceptually  is  the  same  as  nodal  release,  only  the  release 
algorithm  is  different.  To  accomplish  stiffness  reduction,  spring  or  a  combination  of 
spring/gap  elements  are  used  (fig.  6).  The  stiffness  in  the  y  -  direction  is  given  by  the 
spring  constant  of  the  spring  elements.  Crack  growth  is  obtained  by  reducing  the 
spring  constant  of  the  crack  -  tip  element. 

The  nodal  release  method  is  probably  the  most  widely  used  crack  growth  simulation 
technique.  The  crack  is  extended  by  releasing  the  crack  -  tip  node.-  At  the  same  time  a 
reaction  force  is  applied  to  the  released  node  and  then  incrementally  decreased  to  zero. 
The  amount  of  crack  growth,  therefore,  is  restricted  to  the  element  size  per  step. 
Lamain  [31]  stated  that  only  minor  differences  are  observed,  whether  the  reaction 
force  is  applied  proportionally  or  nonproportionally.  During  the  releasing  process, 
the  external  force  can  be  changed  or  held  constant.  It  is  possible  to  use  higher  order 
elements  for  this  method,  but  care  must  be  taken  so  that  crack  face  overlapping  due  to 
the  reaction  force  cannot  occur. 


One  of  the  first  numerical  crack  growth  simulations  was  performed  by  Anderson  [32] 
in  1972,  in  which  he  introduced  the  nodal  release  technique.  For  the  case  of  plane 
stress  and  under  the  assumption  of  a  constant  crack  -  tip  opening  angle  as  the  crack 
growth  criteria,  he  released  the  crack  •  tip  node  and  applied  a  reaction  force  to  this 
node  to  maintain  the  initial  zero  displacement  condition.  Then  he  decreased  the  reaction 
force  in  five  equal  steps,  keeping  the  external  force  constant.  Although  the  assumption 
of  a  constant  CTOA  was  quite  arbitrary  (and  incorrect  for  the  first  few  millimeters  of 
crack  growth),  this  work  can  be  considered  as  the  beginning  of  numerically  stable 
crack  growth  simulation. 

Sorensen  [33]  performed  crack  growth  simulations  for  plane  strain  using  Anderson's 
nodal  release  technique.  He  modeled  crack  extension  for  constant  external  loads 
between  equidistant  nodal  points.  He  discussed  different  possible  fracture  criteria  and 
applied  "a  critical  opening  at  a  small  characteristic  material  distance  from  the  crack  - 
tip’  as  a  criterion  for  stable  crack  growth. 

In  the  1970's,  Shih  et.  al.  performed  extensive  experimental  and  numerical  research 
to  find  valid  crack  growth  criteria  [34],  He  used  the  node  shifting  technique  for  the 
numerical  approach.  His  basic  results  showed  that  the  slope  of  the  J  resistance  curve 
for  A533B  steel  was  constant  for  crack  extension  of  approximately  six  percent  of  the 
remaining  ligament.  Furthermore  he  stated  that  the  "COD  -  based  criteria  appears  to 
be  valid  for  larger  amounts  of  crack  growth’.  The  tearing  modulus  proposed  by  Paris 
et.  al.  [35]  based  on  J  (tearing  modulus:  Tj  -  (E/cr02)(dJ/da))  was  constant  only  for 
a  short  range  of  crack  growth.  An  alternative  approach,  the  tearing  modulus  based  on 
COD  (tearing  modulus:  Td  -  (E/cr02)(d<J/da»  was  considered  to  be  an  ’attractive 
alternative". 

Saka  et.  ai.  [36],  recognizing  the  weakness  of  the  tearing  modulus  concept,  introduced 


in  1983  a  new  tearing  modulus  Tyy  -  (1/RC)  (E/cr02)(dWp/da).  The  parameter  Tw  is 
a  dimensionless  representation  of  dWp/da,  the  incremental  plastic  work  done  in  a 
"circular  region  of  characteristic  radius  Rc*  at  the  growing  crack  -  tip,  where  Tw  is 
directly  related  to  the  amplitude  of  the  singularity  field.  Saka  determined  Rc 
experimentally  to  be  0.28  millimeter  for  A533B  steel  and  performed  a  numerical 
crack  growth  simulation.  The  input  for  the  finite  element  analysis  was  an 
experimental  load  line  displacement  versus  crack  growth  curve.  Saka  compared  Tw 
with  Tj  and  T<j  and  concluded  that  Tw  is  definitely  superior. 

Hoff  [37]  modeled  crack  growth  with  spring  and  gap  elements.  Motivated  by  results 
from  Shih,  he  used  the  J  -  integral  for  only  the  first  four  millimeters  of  crack 
growth.  This  number  was  explicitly  given  by  Kanninen  [18]  as  the  limit  for  J 
controlled  crack  growth  for  A533B  steel.  For  further  crack  extension  Hoff  used  a 
constant  CTOA  value,  verified  by  experimental  data  obtained  from  Shih  [34]. 


Although  much  research  has  been  performed  in  the  numerical  simulation  of  stable 
crack  growth,  up  to  now  an  overall  criterion  that  describes  the  quasi  -  static  extension 
of  a  crack  in  ductile  materials  has  not  be  found.  The  tearing  modulus  concepts  (Tj,  T<j) 
suffer  serious  limitations,  i.e.,  they  show  a  constant  behavior  only  for  a  short  range  of 
crack  growth.  The  same  appears  to  be  true  for  the  new  tearing  modulus  parameter,  Tyy, 

introduced  by  Saka  [36].  In  his  paper  he  performed  a  crack  growth  simulation 
controlled  by  this  parameter  up  to  1.8  millimeter  crack  extension.  His  analysis, 
unfortunately  excluded  the  prediction  of  the  initiation  of  crack  growth.  In  addition,  his 
results  showed  a  deviation  of  about  20  percent  for  Tyy  in  this  crack  growth  interval, 
which  cannot  be  viewed  as  an  improvement  of  over  existing  crack  growth  criteria. 

The  present  work  emphasizes  two  points:  (i)  to  determine  how  well  two  crack  growth 
criteria  (J/CTOD  and  CTOA)  characterize  the  crack  extension  in  mode  1  under  plane 
strain  conditions  and  in  what  interval  range  they  are  applicable,  and  (ii)  to  examine 
the  singular  field  variables  for  crack  growth  initiation  and  subsequent  quasi  -  static 
crack  extension,  since  the  author  believes  that  any  successful  crack  criterion  must  be 
closely  related  to  the  field  variables.  As  input  for  the  finite  element  analysis,  only  an 
experimentally  obtained  load  versus  crack  growth  curve  is  used. 

The  stable  crack  growth  is  simulated  by  using  the  nodal  release  technique.  The  load 
versus  crack  growth  curve  (fig.  7)  is  linearly  discretized  in  a  such  way  that,  when  the 
load  attains  a  certain  value,  the  corresponding  crack  growth  is  modeled  by  releasing  a 
corresponding  number  of  nodes.  The  node  release  is  accomplished  by  replacing  the 
restrained  degree  of  freedom  of  the  crack  -  tip  node  by  a  reaction  force,  which  is  then 
gradually  reduced  to  zero.  After  releasing  the  current  crack  -  tip  node  the  load  again  is 
increased  until  the  requirement  is  satisfied  for  releasing  the  next  node.  The  relation 
used  between  the  applied  load  and  the  crack  growth  for  the  performed  calculations  is 
listed  in  table  1. 
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Since  one  could  argue  that  crack  growth  has  been  observed  to  exhibit  jumping  (pop  - 
in  crack  growth)  behavior  and  that  the  simulation,  therefore,  essentially  is 
linearizing  the  whole  process,  a  second  approach  has  been  performed.  Here  the  node  is 
released  at  a  constant  load  and  the  load  is  increased  with  the  node  restrained  in  the  y  - 
direction.  This  stepwise  or  jumping  simulation  is  shown  in  fig.  8.  The  main 
disadvantage  of  the  latter  simulation  technique  is  the  required  increase  in  CPU  time, 
which  is  nearly  doubled  in  comparison  to  the  first  approach. 

The  crack  tip  opening  displacement  for  both  simulations  are  compared  for  2.25 
millimeters  of  crack  growth.  The  maximum  deviation  occurred  for  the  y  displacement 
at  the  last  node  that  was  released  and  was  always  less  than  three  percent.  For  these 
reasons  the  stepwise  approach  was  not  pursued. 
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The  material  used  for  stable  crack  growth  simulation  is  the  bainitic  pressure  vessel 
grade  A533B  steel.  This  steel  is  representative  of  ductile  materials  and  is  widely  used 
in  pressure  vessel  applications.  The  material  composition  and  the  stress  -  strain 
properties  are  shown  in  table  2  and  table  3.  Ramberg  and  Osgood  proposed  the 
following  constitutive  law  for  simulating  such  material  behavior 


€ 


eq.  5.01 


where  a  and  n  are  material  dependent  constants,  and 

<r0,€0  are  the  yield  stress  and  strain. 


Since  in  the  crack  •  tip  region  the  elastic  strains  are  negligible  in  comparison  to  the 
plastic  strains,  a  simplification  of  eq.  5.01  yields 


€ 


eq.  5.02 


which  is  a  pure  power  law  representation  of  the  stress  -  strain  curve.  Using  a  -  1  and 
n  -  10  the  material  has  been  modeled  with 


€ 


eq.  5.03 


up  to  the  yielding  point  and 


€ 


eq.  5.04 


beyond  yield. 


As  input  for  the  finite  element  program,  the  stress  •  strain  curve  has  to  be 
represented  in  a  multilinear  discretized  form.  This  has  been  achieved  in  16 
discretizations  and  the  actual  input  stress  -  strain  relation  is  shown  in  fig.  9. 

5.2  MODEL  DEFINITION 

A  compact  tension  specimen  is  chosen  to  simulate  mode  1  stable  crack  growth.  This 
specimen  is  chosen  since  it  is  a  standard  type  of  fracture  specimen  that  has  been 
investigated  independently  by  Shih  et  al  [34]  and  Hoff  [37]  and  approximates  plane 
strain  conditions.  The  dimensions  of  the  employed  model  are  shown  in  fig.  10. 

Two  models  are  used  for  the  crack  growth  simulation:  one  coarse  mesh  and  one  fine 
mesh.  The  coarse  mesh  (shown  in  fig.  11a  and  11b)  has  1874  elements  and  1946 
nodes.  The  element  size  in  the  region  of  crack  growth  is  0.25  millimeter.  The  fine 
mesh  shown  in  fig.  12a  to  12c  has  2914  elements  and  3002  nodes.  Here  the  element 
size  in  the  crack  growth  region  is  five  times  smaller  than  in  the  coarse  mesh 
(0.05  mm).  In  the  case  of  the  coarse  mesh,  32  nodes  are  released  which  is  equivalent 
to  eight  millimeters  crack  growth.  For  the  fine  mesh  21  nodes  are  released  to  simulate 
one  millimeter  crack  growth.  For  simplicity,  the  fine  mesh  model  is  created  without 
loading  holes,  however,  previous  work  indicates  that  the  load  is  transfered  to  the  crack 
-  tip  region  reasonably  well  [38,39].  The  boundary  conditions  are  shown  in  fig.  13. 
To  avoid  rigid  body  motion  the  node  ,q  ,  is  restrained  in  x  and  y  -  direction.  All  other 
restrained  nodes  are  restricted  from  moving  only  in  the  y  -  direction. 

5.3  FINITE  ELEMENT  MESH  GENERATION 

An  important  aspect  of  FEM  analysis  is  the  generation  of  the  mesh.  The  two  meshes  used 
for  the  finite  element  analysis,  are  created  with  the  software  package  CAEDS  - 
Graphics  on  an  IBM  5080  terminal.  CAEDS  (  Computer  Aided  Engineering  Design 
System)  is  a  product  of  SDRC  (Structural  Dynamics  Research  Corporation)  employed 
on  a  IBM  4341.  The  strength  of  this  software  package  is  its  flexibility  in  finite 
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element  modeling  and  finite  element  solving.  Since  CAEDS  has  a  direct  interface  with 
CADAM  and  CATIA  (two  engineering  graphic  design  systems),  it  is  possible  to  enter 
part  geometry  directly  from  either  of  those  systems  -  entirely  bypassing  manual 
entry  of  data.  After  doing  that,  one  can  interactively  add  the  load  and  boundary 
conditions  and  submit  the  job  to  the  integrated  finite  element  solver.  The  disadvantage 
of  CAEDS  is  the  restriction  of  its  solver  to  linear  finite  element  analysis.  Other 
problems  encountered  on  CAEDS  are  the  limitation  of  the  number  of  nodes  and  elements 
for  a  successful  analysis  and  the  lack  of  the  support  of  certain  element  types  (like 
plane  strain  elements)  by  the  finite  element  solver.  For  these  reasons,  only  the  mesh 
creations  were  performed  with  CAEDS.  One  of  the  major  challenges  was  to  put  most  of 
the  elements  in  the  vicinity  of  the  growing  crack.  This  has  been  accomplished  by  using 
the  Free  -  Mesh  -  Generator  of  CAEDS.  This  mesh  generator  automatically  creates 
finite  elements  via  the  Triquamesh  algorithm. 


After  a  finite  element  mesh  has  been  created,  the  necessary  mesh  checking  often  is  a 
time  consuming  process.  CAEDS  provides  a  very  powerful  series  of  tools  (the  Model 
Checking  Tools),  for  simplifying  the  model  checking  process.  The  model  checking  tools 
available  in  CAEDS  are: 

*  free  edge  checking, 

-  coincident  node  checking, 

-  interior  element  angle  checking, 

-  distorted  element  checking. 

With  this  checking  series,  not  only  are  modeling  errors  (which  would  make  a  finite 
element  analysis  impossible)  identified,  but  the  elements  are  also  checked  to 
determine,  whether  some  modeling  rules  (aspect  ratio,  element  angles)  are  violated.  A 
violation  of  these  modeling  rules  may  introduce  erroneous  results. 


In  order  to  obtain  improvements  in  the  finite  element  solution  process,  CAEOS  offers 
the  ability  to  optimize  the  bandwith  and  wavefront  of  the  finite  element  mesh.  By 
working  with  free  mesh  generation,  it  is  not  possible  to  number  the  nodes  in  an 
optimum  way.  The  resulting  bandwith  and  wavefront  of  the  final  mesh  thus  becomes 
unacceptably  large  for  meshes  having  2000  •  3000  nodes.  CAEDS  uses  the  Gibbs  - 
Poole  •  Stockmeyer  algorithm  to  optimize  either  the  bandwith  or  wavefront  size  by 
renumbering  the  nodes.  It  is  possible  to  emphasize  the  optimization  either  for  the 
bandwith  or  for  the  wavefront  profile.  After  using  this  optimization  tool  for  the  two 
meshes  created,  the  estimated  CPU  time  reduction  for  performing  the  finite  element 
analysis  was  97  percent.  Without  this  optimization  a  finite  element  analysis  would  not 
have  been  possible. 

After  the  generation  and  optimization  of  the  two  meshes  on  CAEDS,  the  geometry  and  the 
connectivity  of  the  elements  were  transferred  by  a  special  FORTRAN  subroutine  to  the 
finite  element  program  ABAQUS  (40)  for  solution.  This  is  necessary  since  the  CAEDS 
finite  element  solver  does  not  support  plane  strain  elements  (the  state  of  stress  within 
a  compact  tension  specimen  of  this  size  is  assumed  to  be  plane  strain)  and  lacks  the 
ability  to  solve  nonlinear  elastic  or  plastic  problems. 


SJHE  FINITE  ELEMENT  METHOD  FOR  STRUCTURAL  ANALYSIS 

Although  a  considerable  amount  of  work  has  been  performed  to  develop  analytical 
methods  for  solving  problems  of  elasticity  [41,  42,  43]  and  plasticity  [44,  45,  46], 
these  methods  are  usable  only  for  certain  problems  and  analysis  cases.  In  structures  of 
arbitrary  shape  subjected  to  arbitrary  load  conditions,  analytical  methods  often  fail. 
In  engineering  practice,  most  problems  are  too  complicated  be  solved  analytically.  For 
these  cases,  the  finite  element  method  is  a  very  powerful  computational  tool  for 
solving  continuum  mechanics  and  structural  analysis  problems  with  accuracy 
acceptable  to  engineers. 

A  complete  introduction  to  the  finite  element  method  is  far  beyond  the  scope  of  this 
thesis  and,  therefore,  only  a  brief  overview  is  given.  The  interested  reader  is  referred 
to  the  book  of  Zienkiewicz  [47]  which  gives  an  excellent  and  complete  introduction  to 
the  different  approaches  in  finite  element  analysis. 

The  basic  idea  behind  the  finite  element  method  is  to  divide  a  body  into  small 
subvolumes  or  (in  two  dimensions}  a  surface  into  small  subregions.  These  subvolumes 
or  subregions  are  called  elements  and  are  interconnected  at  nodal  points  along  their 
boundaries.  In  the  field  of  solid  materials,  this  method  is  used  to  find  the  stresses  and 
displacements  of  the  structure  being  analyzed. 

In  the  displacement  approach  to  the  finite  element  method  (FEM),  the  displacements  of 
the  nodal  points  are  the  basic  unknown  parameters.  To  approximate  the  displacement 
field  within  each  element,  a  set  of  functions  is  chosen.  These  displacement  functions 
are  called  ’  shape  -  functions'.  The  shape  ■  functions  depend  largely  upon  the  number 
of  nodes  associated  to  each  element  and  the  degrees  of  freedom.  As  a  basic  requirement 
they  must  include  all  possible  rigid  body  displacements  as  well  as  all  appropriate 
strain  states. 

If  these  functions,  in  addition,  satisfy  inter  -  element  compatibility  (which  means  that 
the  highest  derivative  in  the  strain  displacement  relation  must  be  finite)  the 
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displacement  field  will  minimize  the  potential  energy  of  the  system.  Then,  the  finite 
element  solution  represents  an  upper  bound  on  total  potential,  and  the  solution  will 
converge  to  the  true  solution  as  the  mesh  size  is  decreased.  Inter  -  element 
compatibility  may  be  violated  to  produce  reasonable  results,  but  an  upper  bound  on  the 

total  potential  is  no  longer  guaranteed. 

The  strain  within  an  element  can  be  determined  in  terms  of  the  nodal  displacements.  As 
a  final  step,  the  constitutive  properties  of  the  material  will  define  the  state  of  stress 
throughout  the  element  and  on  its  boundaries. 

6.1  THE  FINITE  ELEMENT  FORMULATION 

When  a  body  subjected  to  external  forces  is  in  equilibrium,  the  principle  of  virtual 
work  is  given  by 

<JW|  -  <*WE  eq.  6  01 


where  <5W|  is  the  total  strain  energy 

and  <JWE  is  the  external  work. 


Use  of  the  principle  of  virtual  displacements  gives 


6  W, 


and 


eq.  6.02 


where 


eq.  6.03 


de  is  the  strain  vector  associated  with  virtual  strains, 

o'  is  the  stress  vector, 

b  is  the  body  forces, 

V  is  the  volume, 
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is  the  boundary  where  surface  tractions  are  applied, 
fp  is  the  load, 

<5i i  is  the  virtual  displacement  vector,  and 

T  is  the  surface  traction  vector. 

Substituting  eq.  6.02  and  eq.  6.03  into  eq.  6.01  yields 

J  <5eT  a  dV  -  /  /  <JuTb  dV+J  <5uTTdr  ♦  £  <?uTfp\  -  0  eq.  6.04 
v  l  V  “  r  P-1  -J 

In  a  finite  element  representation  for  solid  materials  the  displacements,  strains  and 
their  virtual  counterparts  can  be  expressed  in  the  following  form 

u  -  N  aj«,  cJu  -  N  da*®  eqs.  6.05a.6.05b 

e-Baj®,  <Jl»Bd£®  eqs.  6.06a,6.06b 

or  in  a  convenient  discretized  form  for  finite  element  applications 


u  -  I  Nj  £j«,  <Ju  «  I  Nj  da,8 

€»SBiaj9,  <5£-IBjdaje 


eqs.  6. 07a, 6. 07b 
eqs.  6. 08a, 6. 08b 


where 


aj®,  daj® 

Nj 

ii 


is  the  i  ih  node, 

are  nodal  displacements  and  their  virtual  counterparts, 
is  the  global  shape  function  for  node  i,  and 
is  the  globai  strain  -  displacement  matrix. 


The  nature  of  N,  and  Bj  is  explained  in  more  detail  in  chapter  6.2.  Using  the  principle 
of  virtual  displacement  and  substituting  eq.  6.05b  and  eq.  606b  into  eq.  6.04  gives 


where  the  transposes  of  the  virtual  strains  and  displacements  are 
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r-'v 


( 


<5eT  -  da°’T  BT 
<JuT  -  <Jae-T  Nt 


eq.  6.10 
eq.  6.11 


In  addition  the  stress  strain  relation  is  defined  as 


cr  -  D  (e  -  ea  )  +<rt 


eq.  6.12 


where 


are  the  initial  strains, 
are  the  initial  stresses,  and 
is  the  constitutive  matrix. 


In  the  case  of  plane  strain  for  linear  elastic  materials,  0  can  be  written  as 


( 1  -  2  u  )  (  1  +  u  ) 


1  -  u  u  o 


u  1  -  u  0 


1-2u 


eq.  6.13 


where  E  is  the  Young's  modulus,  and 

u  is  the  Poisson's  ratio. 


For  elasto  -  plastic  materials  D  is  no  longer  a  matrix  containing  only  elastic  constants. 
Basically  two  new  factors  must  be  introduced: 


1)  a  yield  criterion  (F),  and 

2)  a  hardening  parameter  (k). 
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After  incorporating  F  and  <  into  eq.  6.13  it  can  be  shown  [48]  that  the  constitutive 
matrix,  0,  for  elastic  •  plastic  conditions  is 


(  Da)(  Pa)T 
A+(  Da)T  a. 


eq.  6.14 


where 

a 

represents  the  flow  vector  (a  partial  differentiation  of  the  yield 

criterion  with  respect  to  the  stresses) 

and 

A 

is  a  scalar  term,  obtained  as  the  local  slope  of  the  uniaxial 

stress/  plastic  strain  curve. 

Thus,  by  use  of  eq.  6.14,  eq.  6.13  can  be  rewritten 


£  -Pep(i  -fa)+£ 


eq.  6.15 


for  elasto  •  plastic  conditions.  Substituting  eq.  6.15  into  6.09  yields 
tf jf*’ T|(  J PTPep § d v\  a®  -  jBTpepeadV+jBT£-adV-|NTb  dV 

v  V  /  W  \J  \l 


-I/ldr-I,  NTip].0 


eq.  6.16 


Since  da0'^  is  quite  arbitrary  and  not  necessarily  zero,  the  term  in  the  brackets  must 
be  zero  to  satisfy  eq.  6.16,  or 


(Ja^Bdv) 


a  =  f-  _  f—  +  f ,  +  f—  + 

_5a  JTa  _b  _T  _P' 


eq.  6.17 


where 


Jbtd  €adV 


ifa  J  —  _®P 
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eq.  6.18 
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JiT£'adV 


fb  -  jNTb.dV 
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«r-  J/ldr 


>  -  £1  . 


eq.  6.19 


eq.6.20 


eq.  6.21 


eq.  6.22 


Defining 


I*- !«,-  in ♦  it  *  V  *  if' 

K8-  f  BTD„BdV 
—  J ep  — 


eq.  6.23 


eq.  6.24 


allows  eq.  6.17  to  be  written  as 


K®  a®  -  f®, 


eq.  6.25 


where 


K_®  is  the  stiffness  matrix, 

a®  are  the  nodal  displacements, 

f®  are  the  nodal  forces. 


Equation  6.25  can  be  viewed  as  the  final  representation  of  the  finite  element 
formulation  for  a  solid  material. 


In  ABAQUS,  elements  of  the  type  CPE4  (4  node  bilinear  plane  strain  isoparametric 
elements)  were  chosen  for  the  model.  This  eight  degrees  of  freedom  (d.o.f.)  element 
(two  d.o.f.  for  each  node)  has  the  assumed  displacement  field  [49j 
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or 

u  ■  ai  +  ci2X  +  ci3y  +  a^xy 
v  ■  Q5  +  agx  +  ayy  +  agxy 

where  u  is  displacement  in  x  •  direction,  and 

v  is  displacement  in  y  -  direction 
or  in  a  more  succinct  form 

u-Aan 

where  an  are  the  constant  coefficients  of  the  assumed  polynomial. 
The  displacement  field  also  can  be  written  in  the  form 


y, 


y, 


*ey, 

m 


u_-  N  af 

where  ^1  are  the  shape  functions,  and 

ae  are  the  nodal  displacements. 


eq  6.26 


eq.  6.27 
eq.  6.28 


eq.  6.29 


eq.  6.30 


The  desired  shape  functions  can  be  found  directly  from  Lagrange's  interpolation 
formula,  which  leads  to 


Nt  0  N2  0  N3  0  n4  0 

lvJ 

0  N1  0  N20  N3  0  N4 

V.  V4 


eq.  6.31 


Inserting  the  coordinates  of  the  nodes  gives  the  following  shape  functions  for  the 
rectangular  CPE4  element 


1 

(b  -  x)  (c  -  y) 

eq. 

6.32 

4  b  c 

1 

(b  +  x)  (c  -  y) 

eq. 

6.33 

4  b  c 

1 

(b  +  x)  (c  •  +  y) 

eq. 

6.34 

4  b  c 

1 

(b  *  x)  (c  +  y) 

4  b  c 

eq. 

6.35 

where  fig.  14  identifies  the  parameters  used  . 


The  shape  functions  connect  the  nodal  displacements  with  the  displacement  field. 
Similarly  the  strain  displacement  matrix  B  connects  the  nodal  displacements  with  the 
strain  field.  In  the  case  of  plane  elasticity  B  can  be  written  as 


eq.  6.36 


Incorporating  eqs.  6.32  -  6.35  into  eq.  6.36  gives 


Use  of  eqs.  6.32  to  6.36  permits  the  stiffness  matrix  to  be  evaluated  as 


_K®-j  1  t  dxdy 


-c  -b 


eq.  6.38 


To  evaluate  arbitrary  quadrilateral  CPE4  elements,  an  intrinsic  coordinate  system 
defined  for  each  element  has  to  be  introduced.  In  fig.  15  a  linear  element  is  shown  with 
such  a  natural  coordinate  system.  Axes  6  and  r\  pass  through  the  mid  -  points  of 
opposite  sides  and  the  edges  are  defined  by  6  -  ±  1  and  ri  -  ±  1  regardless  of  how  the 
element  is  oriented  in  the  global  coordinate  system  x,y.  As  a  result  of  this  definition 
node  1  has  the  intrinsic  coordinates  5  -  n  -  *1,  node  2  5-1  and  ri  .  -1  etc.. 

Using  the  discretized  form  of  the  displacements  and  strains  eqs.  6.07  to  6.08,  allows 
the  individual  shape  functions  to  be  written  as 


-  6)  (1  -  n) 


N2--J-(1  +€)  (1  -  n  ) 


(i  *5)  (i  +  n ) 

J  4 


n4-—  (i  -  €)  (i  +  *n ) 


eq.  6.39 


eq.  6.40 


eq.  6.41 


eq.  6.42 


In  general,  u  is  parallel  to  the  x-  axis  and  v  is  parallel  to  the  y  -  axis,  but  they  are  not 
necessarily  parallel  to  5  or  n. 


6.3 
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The  evaluation  of  the  formulated  finite  element  problem  is  done  by  substituting 
eqs.  6.07b  and  6.08b  into  eq.  6.09.  Neglecting  loads  as  well  as  initial  strains  and 
stresses  permits  eq.  6.16  to  be  rewritten  in  the  form 

{  J  BTo-  dV  -  J  N^b  dV  -Llldr  }  -  0  <5.43 

i  y  v  r  j 

Then  the  element  representation  developed  in  chapter  6.2  is  used  to  evaluate  all 
contributions  to  eq.  6.43  separately  for  each  element.  The  displacements  for  each 
element  can  be  obtained  from  eq.  6.07a  as 


u«-INj«aj®. 

For  an  isoparametric  element,  the  x  and 
evaluated  as 


fx»' 

k 

1 

*  ■ 

I 

W 

i  *  1 

eq.6.44 

y  coordinates  within  an  element  can  be 


eq.  6.45 


Then  the  Jacobian  matrix  may  be  evaluated  as 


9x 

ae 

9x 

JtL 

3r| 

viNixe  yiNtve 

y  3N,  xe  y  y ® 
,4*1  ^  1  i  *  1  1 


The  volume  of  each  element  is  given  as 

dV®-t  DetJf  d?dn 


eq.  6.47 
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Following  the  same  approach  as  for  the  displacements,  the  strain  displacement 
relationships  for  each  element  can  be  written 


e_8  - 1  Bj®  aj® 


eq.  6.48 


in  the  case  of  plane  strain 


(&) 


(f-)* 


(&)’ 


eq.  6.49 


and  thickness./,  is  chosen  as  unity.  The  chain  rule  may  be  used  to  evaluate 

3N;  8N;  3N; 

dx  d£  dx  dr]  dx 


3N8  _  3N?  dr] 

By  d%  d y  dy  . 


eq.  6.50 


eq.  6.51 


The  derivatives  (35/3x),(8n/5x),  etc.  can  be  obtained  from  the  inverse  of  the 
Jacobian  matrix,  and  the  stress  -  strain  relationship  for  each  element  can  be  written 


f0  =  Vi«3VIl®-ar 


eq.6.52 


Substituting  eq.  6.52  into  the  first  term  in  the  brackets  of  eq.  6.43,  again  neglecting 
loads  and  initial  strains  and  stresses,  gives  as  the  contribution  from  element  e  to  the 
right  side  of  eq.  6.25 


I  i<ij8  Jj8  -  J  Bj 1  ^  Bj  ^  )  dV  eq.  6.53 

where  Kjj®  is  the  submatrix  of  the  element  stiffness  matrix  JK®.  The  contribution  of 
element  e  to  the  body  force  term  fb  is 

f8-  f  N8,Tb9dV 

V  ” 

and  the  surface  traction  term  is 

<;■  I  dr 

The  actual  integrations  were  evaluated  numerically  in  the  intrinsic  coordinate  system. 
The  most  used  widely  method  (as  in  ABAQUS),  is  the  Gauss  quadrature  method.  The 
submatrix  Kjj®  has  the  form 

i  i 

K®.  =  J  J  Bf  'TCPp  3j8Det  J.8d^  dn  eq.  6.56 

The  nodal  forces  at  node  i  caused  by  the  body  forces  and  surface  tractions  are 


eq.  6.54 


eq.  6.55 


i®  -  'bi8  ♦  *Ti8 


with 


and 


f®  »  f  f  N9'Tb0Det  Jed£ dri 
b ,  i  J  J  — i  —  — 

—  -1-1 

f e  -  J  J  N0’V  Det  J_9d^  dr) 


eq.6.57 


eq.  6.58 


eq.  6.59 


The  strains  and  stresses  are  thus  not  determined  at  the  nodal  points  but  at  the  so  called 
'integration  points'.  In  the  CPE4  elements  these  integration  points  are  marked  in  fig. 
16.  The  Gauss  -  Legendre  method  locates  the  integration  points  in  a  way  that  the 
greatest  accuracy  is  achieved  for  a  given  number  of  them. 


Since  large  deformations  occur  during  the  finite  element  simulation  of  the  stably 
growing  crack,  a  short  introduction  is  given  -about  the  solution  procedure  in  elastic  • 
plastic  finite  element  analysis.  The  form  of  the  stiffness  matrix  given  in  eq.  6.24 
suggests  that  a  straightforward  solution  of  the  finite  element  formulation  may  be 
possible.  This  indeed  is  true  for  the  elastic  case  where  an  explicit  relationship  of  the 
form  of  eq.  6.12  (with  &  -_£(e)  for  nonlinear  elasticity)  is  available.  On  the 
contrary,  such  an  explicit  relationship  is  no  longer  available  for  the  complex  nature  of 
plasticity. 

The  approach  employed  to  solve  the  plasticity  problem  utilizes  the  fact  that  the  matrix 
Dap  (eq.  6.14)  is  known  for  a  certain  stress  value  and  loading  direction,  and  the 
stresses  can  be  integrated  as  shown  in  eq.  6.60  from 


d £  -  DgP  d  €_ 


eq.  6.60 


where  _2ep*  is  known  for  a  certain  stress  value  and  loading  direction. 


A  solution  for  eq.  6.60  can  be  obtained  with  incremental  mathematical  procedures. 
During  the  iteration  process  of  the  elasto  -  plastic  analysis,  the  equilibrium  equation 
(eq.  6.01)  cannot  be  exactly  satisfied,  thus  a  system  of  residual  forces  ¥  will  exist 
such  that 


where  b_  is  the  body  force  vector,  and 

j_  the  external  force  vector. 

Substituting  the  incremental  forms  of  eq.  6.06  and  6.60  into  6.61  for  an  increment  of 
load  gives 

A*  -  K_aAu  -  (Af  +J  h^AddNA  *0  eq.  6.62 

V  V  7 

where  _Ka  is  defined  as  the  element  stiffness  matrix  (eq.  6.23).  With  the  help  of 
incremental  displacements  au,  an  iterative  correction  of 

<Juk  -  [Kf^]-1  A*k  eq.  6.63 

is  calculated  using  the  Newton  -  Raphson  method  where 

<Jjjk  is  used  as  a  corrective  factor. 

After  a  prescribed  number  of  iterations,  the  improved  displacement  is  determined  by 

Auk+1  „  £Uk  <juk  eq.  6.65 

A_u_k+ 1  now  is  resubstituted  into  eq.  6.62  and  the  residual  force  is  calculated.  In 
ABAQUS  the  maximum  residual  force  is  chosen  by  the  user  with  the  parameter 
PTOLyMTOL.  If  the  calculated  residual  force  is  too  high,  the  whole  iteration  process 
must  be  repeated.  The  disadvantage  of  this  procedure  is  that  the  stiffness  matrix  j<a 
must  be  calculated  during  each  iteration.  Therefore  this  method  is  usually  avoided  in 
large  finite  element  codes.  An  alternative  numerical  procedure  is  the  modified  Newton 
Raphson  method,  where  the  stiffness  matrix  is  only  occasionally  recalculated.  The 
initial  stiffness  method  is  such  a  modified  Newton  Raphson  method,  where,  for  the  the 
whole  iteration  process,  the  initial  elastic  stiffness  matrix  is  used. 


The  J  -  integral,  which  is  equal  to  the  strain  energy  release  rate  (for  elastic 
materials),  was  first  introduced  by  Rice  [14]  in  1968.  Under  the  assumption  of  a 
linear  or  nonlinear  elastic  material  free  of  body  forces  and  subjected  to  two 
dimensional  deformation  fields  (i.e.,  plane  strain  or  plane  stress),  the  closed  line 
integral,  J,  around  a  notch  parallel  to  the  x  -  axis  is  path  independent  (see  Fig.  17). 
The  J  -  integral  is  defined  by 


J  -  J  (  Wdy  -  Tj  ds  ) 

eq.  7.01 

r 

where 

€ 

W  »  W(x,y)  -  W(e)  cr.de  . 

eq.  7.02 

o  'J  ,J 

is  the  strain  energy  density  (equivalent  to  the  area  under  the  nonlinear  stress  -  strain 

curve).  Also 

€  "  (€jj) 

eq.  7.03 

is  the  infinitesimal  strain  tensor, 

Tj  ■  (7 ijHj 

eq.  7.04 

is  the  traction  vector  defined  according  to  the  outward  normal  along  r, 

Uj  is  the  displacement  vector,  and 

ds  is  an  element  of  arc  length  along  r. 

The  proof  of  the  path  independence  of  J  is  given  in  (14]  for  a  notch  with  a  finite  radius 
r.  However,  an  application  of  this  integration  formulation  is  also  possible  for  sharp 
cracks,  if  an  arbitrary  small  curve  r  around  the  crack  -  tip  is  assumed,  which 
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reduces  to  zero  in  the  limit.  The  crack  -  tip,  therefore,  can  be  interpreted  as  a 
singularity  of  the  deformation  field.  McMeeking  [50]  determined  the  value  of  J  near 
crack  ■  tip  and  found  out  that  the  evaluation  of  J  in  the  vicinity  of  the  crack  -  tip  is 
not  accurate.  A  practical  limit  on  the  size  of  the  J  -  integral  contour  for  the  mode  1 
compact  tension  specimen  has  been  pointed  out  by  Hoff  [37].  He  suggests  that  J  should 
not  be  calculated  along  contours  closer  than  5 <5  from  the  crack  -  tip,  where  <5  is  the 
crack  tip  opening  displacement. 

For  elastic  -  plastic  calculations  in  the  region  of  small  -  scale  yielding,  the  J  - 
integral  now  is  used  extensively  in  fracture  mechanics,  instead  of  the  stress  • 
intensity  factor  K,  which  is  only  valid  for  linear  elastic  calculations.  In  the  elastic 
case  the  J  •  integral  is  equal  to  the  elastic  potential  energy  release  rate,  G.  By  using 
the  principle  of  virtual  wotk,  it  is  possible  to  derive  the  relation 


eq.  7.05 


where 


as-  1 


as  -(1-  u2) 


for  plane  stress, 

for  plane  strain, 

is  the  Young's  modulus, 

is  the  Poisson’s  ratio, 


and  therefore  (with  J  -  G  in  the  linear  elastic  case) 


J  =  ae 


eq.  7.06 


A  second  method  to  evaluate  the  J  -  integral  is  the  'Virtual  Crack  Extension*  method, 
first  introduced  by  Parks  [51]  for  the  linear  elastic  case  and  later  extended  to  non  - 
linear  material  behavior  [52].  Since  this  method  can  be  implemented  very  easily  in 
FEM  codes,  most  commercial  codes  like  ABAQUS  and  MARC  use  this  technique.  The 
technique  is  based  on  moving  nodes  a  small  distance  around  the  crack  -  tip  by  and 
estimating  the  energy  change.  Since  the  potential  energy  does  not  change  much  with  a 
slightly  different  crack  -  tip  configuration,  this  method  works  very  well.  Problems 
are  encountered,  however,  when  using  collapsed  elements.  In  this  case  only  one  of  the 
several  existing  crack  -  tip  nodes  is  fixed.  The  displaced  crack  configuration  is  so 
different  from  the  original  one  that  the  result  for  the  J  -  integral  would  be  completely 
incorrect.  This  can  be  avoided  by  using  two  rings  of  elements  at  some  distance  away  but 
enclosing  the  crack  -  tip  [31]. 

7.1.2  J  -  1NTEGRALAS  CRITERION  FOR  STABLE  CRACK  -  GROWTH 

As  discussed  earlier,  the  J  -  integral  is  valid  (i.e.,  the  J  -  integral  is  path 
independent)  only  for  elastic  materials  subjected  to  two  dimensional  deformation 
fields. 

Goldman  and  Hutchinson  [53]  showed  that  the  J  -  integral  formulation  may  be  extended 
even  to  elastic  -  plastic  materials  for  cases  of  monotonically  increasing  load  (no 
unloading).  This  implies  that  J  is  strictly  valid  for  analyzing  only  stationary  cracks, 
since  one  of  the  characteristics  of  crack  growth  is  elastic  unloading  and  non  - 
proportional  plastic  deformation  near  the  crack  -  tip.  Nevertheless,  the  J  -  integral  is 
also  used  to  analyze  crack  growth  for  small  amounts  of  crack  extension  [18,  34,  37], 
primarily  because  of  the  lack  of  other  reliable  crack  growth  criteria. 
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Hutchinson  and  Paris  [54]  have  examined  the  necessary  conditions  for  J  controlled 
crack  growth  and  concluded  that  the  most  important  consideration  for  using  J  as  a 
crack  growth  criterion  is  that  nearly  proportional  plastic  deformation  occurs.  In  this 
case,  the  deformation  theory  of  plasticity  and  the  incremental  flow  theory  yield  nearly 
identical  results.  Figure  18  shows  a  typical  J  -  resistance  curve  for  an  intermediate 
strength  steel  under  plane  strain  conditions.  The  dominant  strain  field  as  derived  in  the 
deformation  theory  is 

n 

«lf  kn(-7-r*'?ij<9>  e*707 


where  kn 

is  a  constant. 

r.  9 

are  planar  -  polar  coordinates  centered  at  the 

crack  •  tip 

and  €y 

is  a  function  which  depends  on  n,  the  strain 

hardening  index,  and  whether  plane  stress  or 
plane  strain  is  involved. 

In  Fig.  19  the  crack  -  tip  conditions  are  schematically  shown.  Elastic  unloading  occurs 
only  in  the  direct  vicinity  of  the  crack  extension  zone  (Aa).  However,  it  is  difficult  to 
define  the  size  of  the  zone  where  the  loading  is  nonproportional.  Kanninen  [18] 
suggests  that  this  zone  size  is  of  the  order  of ,/,  which  is  shown  in  Fig.  1 8.  It  should  be 
clear  by  inspection  that  one  condition  for  J  controlled  crack  growth  is  that 


Aa  «  R. 


eq.  7.08 


For  mode  1  the  crack  is  assumed  to  advance  by  an  amount  da  in  the  x  -  direction.  The 
resulting  increment  in  the  strain  field  is 
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For  a  coordinate  system  attached  to  the  crack  tip  using 


3x 


cose 


_a_ 

dr 


sin  e  a 

r  a©  , 


eq.  7.10 


eq.  7.09  becomes 


where 


eq.  7.1 1 


Ve)' 


n 
n  + 


cose€|j(e)+  sine 


a©  €y (e ) 


eq.  7.12 


Inspection  of  eq.  7.1 1  shows  that  the  first  term  corresponds  to  a  proportional  loading 

(for  dJ  >  0)  and  therefore  d€,j  «  €jj.  The  second  term,  however,  is  nonproportional.  It 

is  easy  to  see  that  the  second  term  in  the  bracket  is  of  the  same  order  of  magnitude  as 
the  first  term.  Therefore,  J  controlled  crack  growth  should  be  valid  if  the  proportional 
loading  term  is  much  larger  than  the  nonproportional  loading  term,  or 

dJ.  .da.  eq.7.13 

J  r  . 

By  definition 

J_  =,  4L  _L_  eq.7.14 

I  da  J 

where  /  again  can  be  viewed  as  the  initial  crack  growth  associated  with  the  doubling  of 
J  above  Jc,  fig.  18  [18]. 


then  there  exists  an  annular  region 


l«r «  R  eq. 7.16 

in  which  the  plastic  loading  is  predominantly  proportional  and  the  singularity  field 
(ie.  Eqs.  7.07  and  7.11)  is  dominant. 

By  introducing  a  nondimensional  parameter  defined  by 

U  -  -7-  %  “I-7-'7 

J  03  j 

it  is  possible  to  formulate  a  condition  for  J  controlled  growth  in  a  fully  yielded 
specimen.  Here  b  is  the  uncracked  ligament  and  R  will  be  a  fraction  of  b.  Thus,  finally, 
cj  »  1  can  be  stated  as  requirement  for  J  controlled  growth. 

7.1.3  CALCULATION  QF  THE  J  -  INTEGRAL 

The  integration  of  the  J  -  integral  in  the  CT  specimen  is  performed  along  a  rectangular 
path  which  is  divided  into  six  sections.  Due  to  the  symmetry  of  the  plate  and  loading 
only  the  upper  half  of  the  integral  is  evaluated  and  the  result  is  then  multiplied  by  2. 
The  integration  paths  are  shown  in  fig.  20.  By  use  of  eq.  7.01,  the  J  *  integral  can  be 
separated  into  two  parts 


J  »  Jw  -  J-j- 


eq.  7.18 


The  actual  integration  of  J  was  done  numerically  using  the  average  values  of  the 
stresses,  strains  and  strain  energies  from  the  integration  points  of  the  participating 
elements.  Second  and  fourth  order  finite  difference  operators  were  used  for  evaluating 
3uy/3x. 

To  verify  the  path  independence,  the  J  -  integral  was  calculated  for  4  paths  (  see 
Fig.  21).  The  J  -  integral  values  of  the  4  paths  differ  by  8.2%  which  can  be  viewed  as 
the  deviation  from  the  path  independence  of  J  using  this  numerical  approach  (Table  4). 
For  the  extension  of  the  crack  the  J  -  integral  is  calculated  on  path  2  (  Table  5). 
Figure  22  shows  the  J  -  integral  plotted  against  the  crack  growth.  During  the  first 
millimeter  of  crack  growth  J  exhibits  an  unexpected  behavior.  For  the  first  0.5 
millimeter  of  crack  growth  the  J  value  is  increasing  less  than  for  the  next  0.5 
millimeter.  After  1  millimeter  of  crack  extension  the  J  versus  a  a  curve  shows  the 
expected  behavior.  For  the  first  2.5  millimeters  of  crack  growth  the  J  value  increases 
nearly  linearly;  beyond  that  point  the  slope  is  decreasing.  An  interesting  point  is  that 
the  load  line  -  displacement  versus  crack  growth  (Aa)  curve  shows  a  similar  behavior 
(fig.  23).  By  using  a  linear  least  square  interpolation  for  the  first  2.5  millimeters  of 
crack  growth,  a  slope  of  393.1  MPa  is  obtained.  The  initial  value  of  J  was  89.878 
N/mm  at  a  load  of  5300N. 

Table  6  shows  the  results  of  the  finite  element  analysis  compared  with  the  results 
from  Hoff  [37].  Hoffs  J  versus  A  a  curve  (which  actually  was  his  input  for  the  first 
four  millimeters  crack  growth)  reproduces  the  experimental  data  almost  exactly.  The 
reason  why  the  initial  J  values  differ  so  much  is  easy  explained.  Hoffs  load  ersus  a  a 
curve  (as  a  result  of  his  calculation)  is  significantly  higher  than  the  experimentally 
obtained  curve.  This  however,  was  the  input  of  the  calculation  performed  in  this  paper. 

Therefore,  his  load  at  crack  initiation  was  also  higher  (  -  7000N  )  which  explains  the 
higher  value  of  the  initial  J  -  integral.  The  slopes  of  the  J  versus  Aa  curves  agree 
very  well  with  Hoffs  prediction. 

t 

Using  eq.  7.17  to  determine  whether  J  controlled  crack  growth  is  reasonable  leads  to 
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in  comparison  to  -  150  obtained  by  Hoff.  The  difference  between  these  results  is 
mainly  caused  by  the  different  J  -  values  at  crack  initiation.  Kanninen  (18)  pointed 
out  that  the  question  of  what  is  the  smallest  value  of  u>  for  which  J  controlled  crack 
growth  is  assured  remains  unanswered. 


Although  extensive  research  has  been  done  to  support  the  use  of  the  J  -  integral  as  a 
crack  growth  criterion  for  small  amounts  of  crack  growth,  the  validity  of  such  work  is 
still  doubtful.  The  ui  -  value  (which  should  be  considered  as  the  overall  crack  growth 
criteria)  obtained  in  this  research,  differs  significantly  from  those  reported  in  (18]. 
A  lower  bound  for  u  is  not  known.  In  addition,  the  strain  hardening  exponent  n  and  the 
state  of  stress  has  a  large  influence  on  eq.  7.17.  The  calculation  performed  in  this 
work  and  its  comparison  to  other  results  showed  that  the  slope  (AJ/Aa)  is  most 
reliable  for  the  use  as  a  crack  growth  criteria  for  smaller  amounts  (Aa<  0.036b  )  of 
crack  growth. 


Many  authors  have  employed  fracture  criteria  based  on  the  crack  opening  profile 
geometry.  Examples  include  the  COD  (Crack  Opening  Displacement),  CTOD  (Crack  Tip 
Opening  Displacement)  and  CTOA  (Crack  Tip  Opening  Angle).  Unfortunately,  the 
critical  values  of  these  parameters  are  highly  sensitive  to  their  precise  definition. 
There  is  no  universal  agreement  on  appropriate  definitions  for  COD,  CTOD,  or  CTOA.  in 
fact,  Schwalbe  noted  that,  at  a  recent  conference,  no  less  than  seven  different 
definitions  of  CTOD  were  presented  [55].  The  main  objective  of  using  crack  profile 


m 


geometry  parameters  is  to  describe  the  conditions  at  the  crack  -  tip,  or  to  find 
characteristic  parameters  which  describe  the  crack  -  tip  condition  from  the  beginning 
of  crack  initiation  until  failure,  including  stable  as  well  as  unstable  crack  growth. 


7^LDEFIN1TIQN  OF  THE  CRACK  -  PROFILE  PARAMETERS 

Two  different  crack  profile  parameters  are  studied  in  this  work: 

a)  dt  -  the  'tangent'  definition  of  the  crack  opening  displacement. 

Actually  the  deformed  far  field  crack  front  is  (in  this  case)  extrapolated 
linearly  to  the  original  crack  •  tip  (fig.  24).  This  definition  is  often 
used  with  FEM  calculations  when  the  crack  •  tip  is  not  modeled 
accurately  enough  to  show  that  the  crack  tip  opening  angle  at  crack 
initiation  is  tt  radians.  A  necessary  condition  for  the  use  of  <5X  as  crack 

tip  opening  displacement  is  that  the  deformed  far  field  crack  front  be  a 
straight  line.  In  this  analysis  this  requirement  is  satisfied  for  the  fine 
as  well  as  for  the  coarse  mesh. 

8)  act  -  defined  as  the  crack  tip  opening  angle  (CTOA).  A  commonly  used 
representation  of  act  is 

g-..  tan'  1  2  (Vct  t  l)  eq.  7.25 

h 

vct  is  the  y  -  displacement  of  the  first  node  beyond  the  crack  -  tip,  and 
h  is  the  element  size. 

Since  the  crack  growth  is  simulated  with  uniform  step  size,  the 
determination  of  act  from  this  definition  can  be  quite  accurate. 
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7.22  DISCUSSION  OF  THE  CRACK  PROFILE  PARAMETER 

The  deformed  crack  profiles  for  the  coarse  as  well  as  for  the  fine  mesh  are  shown  in 
fig.  25  to  28.  For  better  display,  the  vertical  scale  is  amplified  in  the  figures.  The 
actual  profiles  are  given  for  both  meshes  in  fig.  29a  to  29q  and  fig.  30a  to  30h.  The 
fine  mesh  predicts  a  blunter  crack  opening  profile  than  the  coarse  mesh,  which  is 
more  consistent  with  experimental  observations. 

L2J3  CRACK.Tlg^g£NMa  DISPLACEMENT  (CTQDi  -  RESULTS 

The  CTOD  versus  crack  growth  (Aa)  diagram  for  the  coarse  mesh  is  shown  in  fig.  31. 
Since  the  J  -  integral  and  the  CTOD  are  similar  concepts,  both  the  J  -  Aa  and  the  CTOD 
•  Aa  curves  have  the  same  shape.  The  slope  of  the  CTOD  •  Aa  curve  increases  for  the 
the  first  millimeter  of  crack  extension,  although  not  to  the  same  extent  as  the  J  -  a  a 
curve.  This  behavior  is  confirmed  by  fig.  32  which  shows  the  CTOD  -  Aa  curve  from 
the  fine  mesh  for  the  first  millimeter  crack  extension.  After  1  millimeter  of  crack 
extension  the  slope  of  the  CTOD  versus  A  a  curve  begins  to  decrease.  The  results 
demonstrate  that  the  CTOD  versus  A  a  curve  is  approximately  linear  for  crack 
extension  up  to  2.5  millimeters.  After  2.5  millimeters  of  crack  growth,  however,  no 
linearity  is  apparent. 

7.2.4  CRACK  TIP  OPENING  ANGLE  (CTO At  -  RESULTS 

The  problem  of  establishing  a  standard  definition  for  the  CTOA  has  been  examined  by 
various  authors  [33,  56],  Rice  [57,  58]  obtained  for  stable  crack  growth  (non  - 
hardening  materials)  a  displacement  distribution  proportional  to  !n(l/r).  Applying 
this  distribution  leads  to  the  conclusion  that  the  crack  tip  opening  angle  is  not  defined 
for  r  »  o  since  d<J/dr  -*  «  as  r  -»  0,  a  result  that  has  been  observed  experimentally 
[18,36]. 

The  CTOA  -  Aa  curve  calculated  for  the  coarse  mesh  is  shown  in  fig.  33.  For  the  first 
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two  millimeters  of  crack  growth,  the  CTOA  -  Aa  curve  shows  a  completely  unstable 
behavior.  It  should  be  clear  that  eq.  7.25  and  the  given  element  size  h  can  not  simulate 
the  CTOA  very  close  to  the  crack  -  tip.  This  definition  should  be  considered  as  a  secant 
approximation.  This  is  especially  true  for  high  strain  hardening  exponents  (n  -  10  in 
our  case).  In  fig.  34  an  attempt  is  made  to  show  the  sensitivity  of  the  CTOA  definition  to 
the  onset  of  crack  growth.  It  can  be  seen  that  the  angles  are  smaller  when  larger 
elements  are  used.  However,  the  inconsistency  in  these  results  was  the  primary  reason 
for  developing  the  fine  mesh  (the  first  CTOA  was  approximately  0.4  radians  in 
comparison  to  a  theoretical  value  of  tt  radians).  The  results  of  the  fine  mesh  (fig.  35) 
show  an  initial  value  of  one  radian  at  the  beginning  of  the  crack  growth  (which  agrees 
with  the  experimentally  observed  crack  blunting  much  better  than  the  coarse  mesh). 
After  the  rapid  decrease  of  the  CTOA  -  a  a  curve  for  the  first  node  release  (0.05  mm 
step  width),  however,  the  fine  mesh  shows  the  same  trend  as  the  coarse. 

Hoff  [37]  also  observed  an  unstable  behavior  in  the  CTOA  -  a  a  curve.  He  stated  as 
reasons  the  ambiguous  definition  of  a  as  well  as  mesh  refinement  errors  due  to  his 
node  release  technique  with  gap  and  spring  elements.  This  is  not  true  for  the 
calculation  presented  in  this  thesis,  since  the  step  -  width  was  constant  during  node 
release.  An  explanation  for  the  present  instability  could  be  that  CTOA  is  not  well 
defined  when  r  approaches  zero.  More  research  is  needed  to  get  further  information 
about  the  nature  of  the  CTOA  for  r  -+  0. 

A  constant  CTOA  -  0.23  radians  is  achieved  after  six  millimeters  of  crack  extension, 
which  is  in  excellent  agreement  of  the  experimental  value  of  0.22  radians  [37]. 


Shih  [59]  developed  a  . t-lationship  between  the  CTOD  and  the  J  -  integral  for  a  static 
crack  that  exploits  the  dominance  of  the  HRR  -  singularity  in  the  crack  -  tip  region. 
He  obtained 
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where  <J45  is  the  crack  tip  opening  displacement  defined  in  fig.  36.  and 

dn  is  a  constant. 

The  constant  dn  depends  mainly  on  the  state  of  stress  and  on  the  strain  hardening 
coefficient  n.  Although  he  used  a  different  definition  for  the  crack  tip  opening 
displacement,  a  comparison  of  the  results  with  the  results  obtained  in  the  fine  mesh 
calculation  is  possible.  The  constant  dn  was  evaluated  from  fig.  37  which  shows  the 

dependence  of  dn  from  n  and  cr<J  E  for  plane  strain. 

Use  of  this  diagram  as  illustrated  gives  dn  *  0.5.  The  calculation  for  J  at  crack 
initiation  yields  102.53  N/mm  in  comparison  to  88.44  N/mm  by  direct  evaluation  of 
J  with  the  line  integral.  Since  it  can  be  expected  that  the  45  degree  definition  of  <5  Shih 
used  in  deriving  eq.  7.28  would  lead  to  a  slightly  lower  value  of  <5,  the  results  can  be 
viewed  as  in  good  agreement. 


8.  INVESTIGATIONS  OF  THE  FIELD  VARIABLES  FOR  STABLE  CRACK  GROWTH 


One  approach  to  representing  the  singular  field  in  the  vicinity  of  the  crack  -  tip  is  to 
describe  it  in  terms  of  the  strain  singularity.  For  stationary  cracks  Rice  and  Rosengren 
[13]  and  Hutchinson  [12]  developed  solutions  for  the  near  crack  -  tip  fields  using  the 
deformation  theory  of  plasticity  and  a  power  hardening  law.  The  stresses,  strains  and 
displacements  are  of  the  form 


°ij-  *0 


(  a€0cr0lnr  ) 


n  +  1 
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eq.  8.01 
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eq.8.02 


eq.  8.03 


where  cr  jj.  €jj  and  Ujj  are  functions  of  e  and  n, 

a  is  the  coefficient  of  the  Ramberg  Osgood  material  description, 
ln  a  constant  given  in  [18], 

cr0,€0  are  the  yield  stress  /  strain. 


A  stably  growing  crack  in  a  ductile  material  causes  large  deformations  in  front  of  and 
elastic  unloading  behind  the  crack  -  tip.  Cracks  opened  by  tensile  mode  1  loadings  are  of 
particular  interest  since  most  fracture  failures  occur  under  mode  1  conditions. 
Unfortunately,  the  mathematical  problems  are  so  complex  that  no  general  analytical 
solutions  are  available. 


The  characteristics  of  the  elastic  -  perfectly  plastic  strain  singularity  of  plane  strain 
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stable  crack  growth  has  been  investigated  by  Rice  (57,58],  Rice  and  Sorensen  [30]  and 
Cherepanov  [60].  In  these  analyses,  the  Prandtl  slip  line  theory  (a  technique  where  the 
stresses  are  determined  by  interpreting  a  slip  line  diagram  according  to  prescribed 
rules)  was  used  to  investigate  the  nature  of  the  elastic  -  plastic  strain  singularity  in  the 
centered  fan  sector  moving  with  a  growing  crack.  The  Prandtl  field  results  for  a 
stationary  crack  have  been  modified  by  introducing  an  elastic  unloading  sector  for  the 
advancing  crack.  It  has  been  found  that  this  sector  is  approximately  between  © ,  -115° 
and  ©2  -  163°  [61].  This  sector  is  shown  in  fig.  38. 

For  a  stationary  crack,  the  elastic  •  perfectly  plastic  asymptotic  strain  singularity  is 
proportional  to  (1/r)  [12,13].  The  results  of  the  investigations  for  plane  strain  stable 
crack  growth  suggest  that  the  nature  of  this  singularity  changes  to  a  weaker  ln(1/r) 
proportionality  [59].  The  reason  for  this  could  be  that  the  crack  is  extending  into 
material  which  has  already  been  deformed  plastically  so  that  complete  refocusing  of  the 
strain  field  ahead  of  the  crack  -  tip  is  prevented  [33]. 

Drugan  et  al.  [62]  constructed  an  exact  solution  for  the  plane  strain  near  -  tip  stress 
field  of  an  advancing  crack  for  nonhardening  materials  by  specializing  Rice's  more 
general  formula.  Sham  [63]  performed  a  finite  element  study  to  verify  this  solution  and 
his  results  agreed  very  well  with  Drugan's  analytical  predictions. 

In  a  more  general  investigation  of  stress  -  strain  fields  for  stably  growing  cracks, 
Amazigo  and  Hutchinson  [64]  examined  a  linear  strain  hardening  material.  Using  J2  flow 
theory  of  plasticity,  they  identified  a  loading  and  an  unloading  zone  near  the  crack  -  tip. 
Nevertheless,  they  did  not  include  a  sector  of  reverse  plasticity  (fig.  38  sector  C)  in  the 
wake  of  the  advancing  crack  found  as  a  result  of  Drugan's  solution.  In  deriving  their 
results,  Amazigo  and  Hutchinson  followed  a  procedure  similar  to  the  HRR  singularity 
approach.  A  nonlinear  stress  rate  function  has  been  generated  to  represent  the  stress  - 
strain  fields  and  the  order  of  the  elastic  unloading  zone  has  been  determined  dependent  on 
the  slope  of  the  linearly  simulated  plastic  pad  of  the  stress  -  strain  relation. 


As  mentioned  above,  there  is  no  exact  solution  available  to  describe  the  field  /ariables  of 
a  stably  growing  crack  for  the  general  case  of  a  power  law  hardening  material.  In  the 
discussion  of  the  field  quantities  from  the  present  finite  element  solution,  the  emphasis 
is  placed  on  the  transition  between  the  different  parameter  fields  displayed  by  the 
curves.  The  definition  of  the  representation  of  the  curves  are  given  in  fig.  39.  It  is  clear 
that  for  pure  mode  1  loading  the  line  ©  -  0°  best  shows  the  phenomena  described  above 
and  thus  the  various  quantities  are  shown  along  this  line. 

The  strong  HRR  field  strain  singularity  for  the  stationary  crack  is  quite  evident  on  the 
y  -  direction  strain  curve  along  the  line  ©  -  0®  of  the  fine  mesh  (fig.  40).  It  is  self 
evident  that  the  coarse  mesh  is  not  able  to  simulate  this,  nearly  1/r  singularity  for  the 
stationary  crack.  As  the  crack  grows  larger  and  larger  the  predominance  of  the  HRR  field 
near  the  crack  •  tip  decreases.  Nevertheless,  the  characteristics  of  the  y  -  direction 
strain  curves  change  very  slowly  (fig.  41  to  43).  Only  after  five  millimeters  of  crack 
extension  is  the  strain  field  singularity  observed  to  be  significantly  weaker  than  at  the 
onset  of  crack  growth  (fig.  43).  One  possible  explanation  for  this  strain  curve  behavior 
could  be  the  smoothing  effect  of  the  elastic  strain  component  for  the  first  increment  of 
crack  extension. 

The  y  -  direction  stresses  relative  to  the  crack  -  tip  along  the  line  ©  -  0®  exhibit  little 
change  for  the  first  millimeter  of  crack  growth,  for  the  fine  as  well  as  the  coarse  mesh 
(fig.  44  to  47).  This  is  not  very  surprising  since  the  known  theories  assume  only 
slightly  different  singularities  for  the  near  -  tip  field  of  the  advancing  crack  and  the 
stationary  crack.  For  larger  amount  of  crack  extension  the  singularity  tends  to  become 
stronger.  The  points  marked  in  fig.  45  could  be  interpreted  as  the  transition  between  an 
intermediate  zone,  where  the  plastic  strain  is  comparable  in  magnitude  with  the  elastic 
strain  and  the  K  field.  The  transition  between  the  intermediate  zone  and  the  HRR  field  is 
not  resolvable.  For  a  crack  extension  of  more  than  one  millimeter  the  influence  of  the  K 
field  is  diminished  and  the  transition  points  are  no  longer  identifiable. 


The  von  Mises  stress  tor  the  line  ©  -  0°  clearly  reveals  the  transition  between  the 
intermediate  zone  and  the  K  field  for  both  meshes  (fig.  48  to  51).  The  transition  points, 
which  were  hardly  detectable  for  the  stress  in  y  -  direction  (fig.  45),  are  now  very 
evident.  Since  the  applied  external  load  increases  as  the  crack  advances,  the  transition 
points  change  toward  larger  r  values.  The  crack  growth  versus  transition  point  position 
is  depicted  in  fig.  52.  At  the  beginning  of  crack  extension  the  transition  points  change 
very  rapidly,  while  the  change  decreases  as  the  crack  grows  larger.  For  the  growing 
crack,  a  resolution  of  the  crack  -  tip  singularity  from  the  HRR  field;  as  well  as  the  HRR 
field  from  the  intermediate  zone  is  not  possible  even  with  the  fine  mesh  (fig.  51).  For 
closer  examination,  a  log/log  representation  of  the  von  Mises  curve  is  plotted  in  fig.  53 
and  54.  Indeed,  the  changes  in  slope  circled  in  fig.  53  could  be  interpreted  as  the 
transition  between  the  HRR  field  and  the  intermediate  zone. 

The  elastic  strain  energy  density  (  We| )  is  shown  for  the  early  stages  of  crack  extension 
in  fig.  55  for  the  fine  mesh  along  the  line  ©  *  0°  &  1 80°  (negative  values  of  r  represent 
©  -  180°).  After  the  onset  of  crack  growth,  the  elastic  strain  energy  density  is  seen  to 
decrease  dramatically.  For  the  next  steps,  the  magnitude  of  We|  was  nearly  constant, 
although  the  applied  external  load  increased  significantly  during  crack  extension.  This 
behavior  characterizes  in  an  excellent  way  the  ductile  material  behavior,  namely,  if  the 
the  load  were  to  be  held  constant  instead  of  being  increased,  the  crack  would  arrest.  Due 
to  the  strong  gradient  of  We)  for  the  first  step,  this  behavior  could  not  be  resolved  with 

the  coarse  mesh.  Figures  56  and  57  depict  the  We)  characteristics  for  subsequent  crack 
extension  with  the  coarse  mesh,  and  show  that  the  magnitude  of  the  maximum  value  of 
WQ|  decreases  for  the  growing  crack.  The  maximum  values  of  Wa)  plotted  in  fig.  55  to 
57  should  only  be  interpreted  qualitatively  due  to  the  numerical  inaccuracy  associated 
with  the  elements  located  directly  at  the  crack  -  tip. 

The  plastic  strain  energy  density  (  Wpi  )  for  the  advancing  crack  is  plotted  for  the  line 
9  -  0°  &  180°  in  fig.  58  to  60.  The  curves  show  that  the  maximum  value  of  Wp,  is 
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located  behind  the  crack  -  tip.  It  is  more  important,  however,  that  the  energy 
dissipation  is  sharply  bordered  and,  for  larger  amounts  of  crack  growth,  is  nearly 
constant  over  a  certain  length  that  is  extending  with  the  advancing  crack,  a  result  which 
could  possibly  be  employed  in  a  local  crack  characterization.  As  already  discussed  in 
chapter  2,  Saka  [35]  performed  detailed  research  into  the  feasibility  of  the  plastic 
dissipation  as  a  local  crack  growth  criterion.  He  stated  that  the  intense  region  of  the 
plastic  dissipation  is  circular  with  its  origin  at  the  tip.  In  addition,  he  determined  the 
characteristic  radius  of  this  circle  to  be  0.28  millimeter  for  approximately  two 
millimeters  of  crack  extension.  In  contrast,  the  results  of  the  present  analysis  indicate 
that  the  intense  region  of  Wpj  is  located  behind  the  crack  -  tip.  The  radius  of  this  intense 

region  (assumed  to  have  a  circular  shape)  has  been  determined  to  be  approximately  one 
millimeter  for  two  millimeters  crack  extension  and  is  clearly  a  function  of  the  crack 
extension. 

For  the  estimation  of  the  plastic  zone  size,  the  von.  Mises  stress  is  illustrated  for  the 
advancing  crack  in  fig.  61a  to  61m.  In  the  present  analyses  the  yield  stress  is  382  MPa. 
At  the  onset  of  crack  extension  the  plastic  zone  exhibits  the  typical  butterfly  shape,  and, 
as  the  crack  grows,  the  plastic  zone  becomes  more  characterized  by  a  bending  behavior 
towards  the  uncracked  ligament.  After  Q.25  millimeter  of  crack  growth,  plastic  hinging 
occurs  at  the  end  of  the  uncracked  ligament  due  to  the  moment  caused  by  the  applied 
external  load.  At  step  4  (Aa  -  0.75  mm,  P  *  7650  N)  the  plastic  zone  of  the  crack  joins 
the  plastic  hinging  region.  With  increasing  crack  growth,  the  crack  -'  tip  moves  from  the 
middle  of  the  highest  von  Mises  stress  contour  to  its  left  border  and  it  seems  that  for 
further  crack  growth  the  crack  -  tip  would  restrain  the  highesl  von  Mises  stress 
contour.  One  interesting  region  of  the  compact  tension  specimen  is  right  in  front  of  the 
crack  -  tip.  In  this  region  the  magnitude  of  the  von  Mises  stress  is  decreasing  towards 
the  line  ©  -  0°.  A  reasonable  explanation  for  this  could  be  the  superposition  of  the  stress 
parallel  to  the  external  load  with  the  stress  due  to  the  resulting  moment. 
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The  theory  of  stable  crack  growth  in  ductile  materials  suggests  that  the  region  ahead  a 
crack  •  tip  should  be  divided  into  four  regions: 

1  )  The  crack  •  tip  singularity  for  advancing  cracks, 

2  )  the  HRR  field. 

3 )  an  intermediate  zone  where  the  plastic  strain  is  comparable  in 
magnitude  with  the  elastic  strain,  and 

3  )  the  K  field. 

In  the  present  finite  element  study,  only  the  transitions  between  the  HRR  field  • 
Intermediate  zone  and  intermediate  zone  -  K  field  could  be  identified.  Neither  the  fine 
nor  the  coarse  meshes  are  able  to  resolve  the  crack  -  tip  singularity  of  the  advancing 
crack  from  the  HRR  field.  Although  one  could  argue  that  the  finite  element  formulation  is 
smoothing  the  results  and,  therefore,  such  a  separation  may  not  be  observable  with  this 
method,  the  nature  of  the  obtained  stress  curves  clearly  suggest  that  the  crack  •  tip 
singularity  is  superposed  onto  the  HRR  field  rather  than  appearing  separately. 


9.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 


In  the  present  work,  slow  stable  crack  growth  in  ductile  material  has  been  simulated 
numerically  with  a  widely  used  commercial  finite  element  code  (ABACUS).  An 
experimentally  obtained  load  versus  crack  growth  relation  was  used  as  input  and 
J/CTOD  and  CTOA  crack  growth  criteria  in  ductile  fracture  mechanics  were 
investigated.  In  addition  the  field  parameters  ahead  of  the  tip  of  the  growing  crack  were 
investigated. 

Existing  theories  about  the  asymptotic  fields  ahead  of  a  crack  -  tip  for  an  advancing 
crack  indicate  that  there  should  exist  four  regions.  The  results  of  the  finite  element 
analysis  show  the  near  crack  •  tip  field  of  the  advancing  crack  and  the  HRR  field 
characterizing  the  stationary  crack  are  superposed  on  each  other  and  do  not  appear 
independently.  The  only  transitions  which  are  resolvable  are  the  transitions  between 
the  HRR  field  -  intermediate  zone  (the  zone  where  the  plastic  strain  is  comparable  in 
magnitude  with  the  elastic  strain)  and  the  intermediate  zone  -  K  field.  The  strain  field 
tends  to  a  significantly  weaker  singularity  as  the  crack  grows  larger,  whereas  the 
stress  field  remains  nearly  unchanged  even  for  large  amounts  of  crack  growth.  In 
addition,  the  plastic  dissipation  energy  field  attains  its  maximum  value  behind  the 
crack  -  tip  and,  for  larger  amounts  of  crack  growth,  the  plastic  energy  dissipation 
possesses  a  nearly  constant  value  over  a  certain  fixed  distance. 

The  results  of  the  finite  element  analyses  performed  in  this  work  show  that  the  J  - 
integral  /  CTOD  concepts  do  not  appear  promising  as  crack  growth  criteria.  This  is  not 
very  surprising  since  the  J  -  integral  concept  is  strictly  valid  only  for  stationary 
cracks,  although  the  slope  of  the  J  resistance  curve  (  aJ/  Aa)  does  appear  to  be  useful 
as  a  crack  growth  controlling  parameter  for  limited  amount  of  crack  growth  (Aa  < 
0.036b).  This  indeed,  is  a  very  important  result  because  other  crack  growth  criteria 
available  are  not  able  to  simulate  this  region  of  crack  growth  very  well.  One  of  these 
criteria  is  the  crack  tip  opening  angle  (CTOA),  in  which  a  constant  CTOA  is  achieved 
after  six  millimeters  of  crack  extension.  Evaluation  of  the  finite  element  results  for 
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the  advancing  crack  indicate  that  the  use  of  the  CTOA,  as  a  crack  growth  criteria,  is  not 
very  accurate  for  small  amounts  of  crack  extension.  In  addition,  any  definition  of  the 
CTOA  fails  for  r  approaching  the  limit  r  -*  0.  Obviously  there  exists  a  gap  between  the 
region  where  the  J  -  integral  is  path  independent  {and  the  slope  of  the  J  -  resistance 
curve  characterizes  the  crack  growth)  and  the  region  where  the  CTOA  is  constant  and 
therefore  applicable  as  a  crack  growth  criterion.  For  this  reason,  a  combination  of 
(AJ/  Aa)  and  the  CTOA  as  the  crack  growth  criteria  proposed  by  Kanninen  and  Popelar 
[18]  and  recently  applied  by  Hoff  [37]  is  at  best,  an  approximation  where  the  extent 
of  the  error  remains  to  be  determined.  To  make  numerical  crack  growth  simulation 
techniques  applicable  for  practical  problems  future  research  should  be  focused  on 
three  points: 


1  )  The  amount  of  crack  growth  for  which  the  J  *  integral  is  nearly  path 

independent  needs  to  be  known  for  a  much  wider  range  of  materials  and 
geometries.  The  w  approach  of  Hutchinson  and  Paris  does  not  seem  promising. 
In  addition,  Kanninen  [18]  pointed  out  that  the  question  of  determining  the 
smallest  value  of  w  to  assure  J  controlled  crack  growth  remains  unanswered. 

2  )  An  unambiguous  definition  of  the  CTOA  needs  to  be  established,  and  the  amount 

of  crack  growth  when  the  CTOA  begins  to  be  constant  needs  to  be  better 
understood.  Since  the  CTOA  appears  to  be  material  dependent  [33],  more 
experimental  research  needs  to  be  done  to  answer  these  questions. 

3  )  An  additional  crack  growth  criterion  needs  to  be  developed  for  the  region 

between  J  -  and  CTOA  controlled  crack  growth.  This  new  criterion  could  even 
replace  one  or  both  of  the  crack  growth  criteria  mentioned  above.  Two  new 
energy  based  crack  growth  criteria  may  be  able  to  satisfy  the  requirement(s) 
mentioned  above: 

(i)  the  plastic  dissipation  energy  criteria  introduced  by  Saka  [36].  The 
idea  behind  this  criterion  is  that  the  plastic  dissipation  energy 


determined  in  an  intense  strain  region  with  a  characteristic  radius  Rc  and 

written  in  a  dimensionless  representation  causes  crack  growth  when  a 
critical  value  is  exceeded.  As  discussed  earlier,  the  results  given  in  Ref.  [36] 
differ  from  the  analysis  presented  in  this  work.  Nevertheless  fig.  62  shows 
how  the  characteristic  radius  Rc  for  the  strain  intense  region  could  be  defined 
from  the  point  of  view  of  this  work.  An  interesting  fact  is  that  for  larger 
amounts  of  crack  growth  the  characteristic  diameter  (2RC)  appears  to  be 
approximately  of  the  order  of  the  crack  extension  Aa. 

(ii)  the  strain  energy  density  criterion  proposed  by  Sih  [65].  The  basic 
hypothesis  behind  this  criterion  is  that  the  maximum  yielding  is  assumed 
to  coincide  with  maximum  strain  energy  density  and  the  fracture  initiation 
with  minimum  strain  energy  density.  Failure  occurs  now  if  either  the 
maximum  strain  energy  density  or  the  minimum  strain  energy  density  exceed 
a  critical  value.  Since  the  strain  energy  density  is  easy  to  measure  (crea 
under  the  stress/strain  curve)  this  new  approach  looks  promising. 


[1]  C.E.  Inglis,  Trans.  Roy.  Inst.  Naval  Architects  60,  219,  (1913). 

(2]  H.  Neuber,  "Theory  of  Notch  Stresses",  Edwards,  Michigan,  (1946). 

“Kerbspannungslehre.  Grundlagen  fuer  genaue  Spannungsrechnung". 
Springer  Verlag,  Berlin ,(1937). 

[31  A.  A.  Griffith,  Phil.  T rans.  Roy.  Soc.  London,  Ser.  A221 ,1 63,(1 921 ). 

[4]  A.  A.  Griffith,  In  Proceedings  of  the  1  st  International  Congress  for  Applied 

Mechanics,  Delft,  p.55.  (1924). 

[5]  I.  N.  Sneddon,  "The  Distribution  of  Stress  in  the  Neighborhood  of  a  Crack  in  an 

Elastic  Solid",  Proceedings  of  the  Royal  Society,  A,  Voi.  1 87,  pp229  - 
238,(1946). 

[6]  H.  M.  Westergaard,  "Bearing  Pressure  and  Cracks",  Transactions  of  the  AMSE, 

Journal  of  Applied  Mechanics,  Vol.  66,  pp.  A49  -  A53,  (1939). 

[7]  G.  R.  Irwin,  J.  Appl.  Mechanics,  24, 361,  (1957). 

[8]  G.  R.  Irwin,  "Handbuch  der  Physik",  Vol.  79,  Springer- Verlag,  Berlin, 

pp.  551  -  590,  (1958). 

[91  E.  Orowan,  "Energy  Criteria  of  Fracture",  Welding  Research  Supplement, 

Vol.  20,  p.  1575,  (1955). 

[1 01  G.  R.  Irwin, "  Relation  of  stresses  near  a  crack  to  the  crack  extension  force", 

Proc.  9th  Int'l  Congress  of  Applied  Mech.,  Uni.  of  Brussels,  Vol.  VIII, 
Chap.  XII,  pp.  245  -  251 ,  (1957). 

[ill  A.  A.  Wells,  "Application  of  Fracture  Mechanics  at  Beyond  General  Yielding", 
British  Welding  Journal,  10,  pp.  563  -  570,  (1963). 

[1 2J  J.  W.  Hutchinson,  "Singular  Behavior  at  the  End  of  a  Tensile  Crack  in  a 

Hardening  Material",  Journal  of  the  Mechanics  and  Physics  of  Solids, 

16,  pp.  13-31,  (1968). 

[1 3)  J.  R.  Rice  and  G.  F.  Rosengren,  "Plane  Strain  Deformation  Near  a  Crack  Tip  in  a 

Power-Law  Hardening  Material",  Journal  of  the  Mechanics  and  Physics 
of  Solids,  35,  pp.  379  -  386,  (1968). 

* 

[14]  J.  R.  Rice,  "A  Path  Independent  Integral  and  the  Approximate 'Analysis  of  Strain 

Concentrations  by  Notches  and  Cracks",  Journal  of  Applied  Mechanics, 
35,  pp.  379-386,  (1968). 


[15]  J.  L  Sanders,  "On  the  Griffith-Irwin  Fracture  Theory",  Journal  of  Applied 
Mechanics,  27,  pp.  352  -  353,  (1960). 


[1 6]  J.  D.  Eshelby,  "Energy  Relations  on  the  Energy-Momentum  Tensor  in 

O'Minuum  Mechanics",  In  Elastic  Behavior  of  Solids,  M.F.  Kanninen  et. 
al.  (ed.),  McGraw-Hill,  New  York,  pp.  77  -  1 15,  (1969). 

[1 7]  G,  P.  Cherepanov, "  On  Crack  Propagation  in  Solids",  International  Journal  of 

Solids  and  Structures,  5,  pp.  863  -  871,  (1969). 

[18]  M.  F.  Kanninen  and  C.  H.  Popelar,  "Advanced  Fracture  Mechanics’,  Oxford 

Engineering  Science  Series,  15,  (1985). 

[19]  J.  F.  Knott,  "Microscopic  Aspects  of  Crack  Extension",  In  Advances  in  Elasto 

Plasic  Fracture  Mechanics,  L.  H.  Larsson  (ed.),  Applied  Science 
Publishers  LTD,  London,  p.  21,  (1979). 

[20]  C.  Zener,  "Micromechanism  of  Fracture",  ASM  40,  p.  3131,  (1948). 

[21]  A.  N.  Stroh,  "The  Formation  of  Cracks  as  a  Result  of  Plastic  Row",  Proc.  Roy. 

Soc.  A  223,  pp.  404  -  414,  (1954). 

[22]  W.  Dahl  and  D.  Dormagen,  "Micromechanisms  of  Crack  Initiation  and  Crack 

Propagation",  In  Porceedings  of  the  4th  Advanced  Seminar  on  Fracture 
Mechanics,  p203,  (1983). 

[23]  A.  H.  Cottrell,  "Theory  of  Brittle  Fracture  in  Steel  and  Similar  Materials", 

Trans.  AIME  212,  p.  192,  (1958). 

[24]  J.  F.  Knott, "  Micromechanisms  of  Fibrous  Crack  Extension  in  Engineering 

Alloys',  Met.Sci.  14,  p.  375,  (1980). 

[25]  F.  A.  McClintock,  In  Physics  of  Strength  and  Plasticity,  A.S.Argon  (ed.), 

M.l.T.  Press,  Cambridge,  Ma,  p.  307,  (1964). 

[26]  J.  R.  Rice  and  Tracey,  Journal  of  the  Mechanics  and  Physics  of  Solids,  Vol.  17, 

p.  201,(1969). 

[27]  J.  R.  Rice  and  Johnson,  In  ’nelastic  Behavior  of  Solids,  M.A. Kanninen  et.  al. 

(ed.),  McGraw-Hill,  New  York,  p.641,  (1970). 

[28]  R.  O,  Richie,  "Slow  Crack  Growth:  Macroscopic  and  Microscopic  Aspects", 

R.B.Tait  and  G.G. Garret  (ed.),  Pergamon  Press  Inc.,  New  York,  p.  100, 
(1984). 

[29]  G.  Green  and  J.  F.  Knott,  Journal  of  the  Mechanics  and  Physics  of  Solids, 

Vol.  23,  p.167,  (1975). 


[30]  J.  R.  Rice  and  E.  P.  Sorensen,  "Continuing  Crack  Tip  Deformation  and  Fracture 

for  Plane-Strain  Crack  Growth  in  Elastic-Plastic  Solids",  Journal  of 
the  Mechanics  and  Physics  of  Solids.  Vol.  26,  pp.  163  -  186,  (1978). 

[31]  L.  G.  Lamain,  "Numerical  Analysis  in  EPFM",  In  Elastic- Plastic  Fracture 

Mechanics,  LH.Larsson  (ed.),  D.  Reidel  Publishing  Company, 
pp.  227  -  261,(1983). 

[32]  H.  Anderson,  "A  Finite  Element  Representation  of  Stable  Crack  Growth", 

Journal  of  the  Mechanics  and  Physics  of  Solids,  Vol.  21 .  pp.  337  -  356, 
(1972). 

[33]  E.  P.  Sorenson,  "A  Numerical  Investigation  of  Plane  Strain  Stable  Yielding 

Condition",  1 7  ASTM  STP  668,  pp.  1 51  -  1 74,  (1979). 

[34]  Shih,  de  Lorenzi  and  Andrews,  "Studies  on  Crack  at  Initiation  and  Stable  Crack 

Growth".  ASTM  STP  668,  pp.  65-  120,  (1979). 

[35]  P.  C.  Paris  et.  al.  "Theory  of  Instability  of  the  Tearing  Mode  of  Elastic-Plastic 

Crack  Growth",  ASTM  STP  668,  pp.  5  -  36,  (1979). 

[36]  M.  Saka  et.  al.,  "A  Criterion  Based  on  Crack  Tip  Energy  Dissipation  in 

Plane-Strain  Crack  Growth  under  Large  Scale  Yielding",  ASTM 
STP  803,  pp.  I  130  - 1 158,  (1983). 

[37]  R.  Hoff,  "A  new  Finite  Element  Technique  for  Modeling  Stable  Crack  Growth", 

Engineering  Fracture  Mechanics,  Vol.  23,  No.  1 ,  pp.  1 05  - 1 1 8, 

(1986). 

[38]  P.  D.  Hilton  and  G.  C.  Sih,  In  Mechanics  of  Fracture,  Noordhoff,  Leyden, 

Netherlands,  Vol.  1 ,  pp.  426  -  476,  (1 973). 

[39]  J.C.  Newman, Jr.,  "Stress  Analysis  of  the  Compact  Tension  Specimen  Including 

the  Effects  of  Pin  Loading",  Fracture  Analysis,  ASTM  STP  560,  ppi05  - 
121,(1974). 

[40]  ABAQUS  Version  4/5/1 75,  Hibbitt,  Karlsson  and  Sorenson,  Inc.,  (1 985). 

[41]  N.  I.  Muskhelishvili,  "Mathematical  Theory  of  Elasticity",  Noordhoff 

International  Publishing,  Leyden,  (1977). 

[42]  I.  S.  Sokolnikoff,  "Mathematical  Theory  of  Elasticity".  Robert  E.  Krieger 

Publishing  Company,  Florida,  (1983). 


[43]  S.  P.  Timoshenko,  "Theory  of  Elasticity",  McGraw-Hill  Book  Company ,(1970). 


[44]  P.  Ludwik,  "Elemente  der  Technologischen  Mechanik",  Springer  Verlag,  Berlin 
(1909). 


[45]  L.  Prandtl,  "Zeitschrift  fur  angewandte  Mathematik  und  Mechanik",  Voi.  8, 

p.85,  (1928). 

[46]  Bauschinger,  Der  Zivilingenieur,  pp.  27  -  289,  (1881). 

[47]  0.  C.  Zienkiewicz,  "The  Finite  Element  Method",  McGraw-Hill,  UK,  3rd  edition, 

(1977). 

[48]  D.  R.  J.  Owen  and  A.  J.  Fawkes,  "Engineering  Fracture  Mechanics",  Pineridge 

Press  Ltd.,  Swansea,  U.K.,  (1983). 

[49]  R.  D.  Cook,  "Concepts  and  Applications  of  Finite  Element  Analysis",  John  Wiley 

and  Sons,  New  York,  (1981). 

[50]  R.  M.  McMeeking,  "Path  Dependence  of  the  J-Intergral  and  the  Role  of  J  as  a 

Parameter  Characterizing  the  Near-Tip  Field",  ASTM  STP  631 , 
pp.  28-41.(1977). 

[51]  D.  M.  Parks,  "A  Stiffness-Derivative  Finite  Element  Technique  for 

Determination  of  Elastic  Crack  Tip  Stress  Intensity  Factors',  Int.  J.  of 
Frac.,  Vol.  10,  pp.  487  -  502,  (1974). 

[52]  D.  M.  Parks,  "The  Virtual  Crack  Extension  Method  for  Non-Unear  Material 

Behavior",  Computer  Methods  in  Mechanics  and  Engineering,  Vol.  12. 
pp.  353  -  364,(1977). 

[53]  N.  L.  Goldman  and  J.  W.  Hutchinson,  Int.  J.  Solids  Structures,  pp.  575  -  591 , 

(1975). 

[54]  J.  W.  Hutchinson  and  P.  C.  Paris,  "Stability  Analysis  of  J-Controlled  Crack 

Growth",  In  Elastic  Plastic  Fracture,  ASTM  STP  668,  pp.  37  -  64, 

(1 979). 

[55]  K.  H.  Schwalbe,  "The  Crack  Tip  Opening  Displacement  in  Elastic  -  Plastic 

Fracture  Mechanics',  Springer-Verlag,  Berlin,  (1986). 

[56]  R.  M.  McMeeking,  "Journal  of  the  Mechanics  and  Physics  of  Solids",  Vol.  25, 

pp.  357  -  381,(1977). 

[57]  J.  R.  Rice,  In  "Fracture:  An  Advanced  Treatise",  H.  Uebowitz,  (ed.)  Vol.  2, 

Academic  Press,  New  York.  pp.  191  -311,  (1968). 


[58j  J-  R-  Rice,  In  "Mechanics  and  Mechanisms  of  Crack  Growth",  M.  J.  May  (ed.), 
British  Steel  Corporation  Physical  Metallurgy  Centre  Publication, 
pp.  14  -39,  (1975). 

[59]  C.  F.  Shih,  "Relationships  between  the  J-Integral  and  the  Crack  Opening 
Displacement  for  Stationary  and  Extending  Cracks",  J.  Mech.  Phys. 
Solids,  Vol.  29.  No.  4.  pp.  305  -  326,  (1981). 

[60j  G.  P.  Cherepanov,  "Mechanics  of  Brittle  Fracture  (in  Russian),  Nauka, 

Moscow,  1974,  English  translation,  Eds.  McGraw-Hill,  New  York, 
pp.  284-289,(1979). 

[61  ]  J.  R.  Rice,  W.  J.  Dragan  and  T.  L.  Sham,  "Elastic- Plastic  Analysis  of  Growing 
Cracks",  In  ASTM  STP  700,  pp.  189  -  221,  (1980). 

[62]  W.  J.  Dragan,  J.  R.  Rice,  and  T.  L  Sham,  J.  Mech.  Phys.  Solids.  Vol.  30, 

p.  447,  (1982)  and  erratum,  Vol.  31,  p.191,  (1983). 

[63]  T.  L  Sham,  “A  Finite-Element  Study  of  the  Asymptotic  Near-Tip  Fields  for 

Model  Plane-Strain  Cracks  Growing  Stably",  In  Elastic  -  Plastic 
Fracture,  ASTM  STP  803,  Vol.  1,  pp.  I  52  - 1  79,  (1981). 

[64]  J.  C.  Amazigo  and  J.  W.  Hutchinson,  "Crack  Tip  Fields  in  Steady  Crack-Growth 

with  linear  Strain-Hardening",  J.  Mech.  Phys.  Solids,  Vol.  25, 
pp.  81  -97,(1977). 

[65]  G.  C.  Sih  and  E.  Madenci,  "Fracture  Initiation  under  Gross  Yielding:  Strain 

Energy  Density  Criterion",  Engineering  Fracture  Mechanics,  Vol.  18, 

No.  3,  pp.  667-677,(1983). 

[66]  W.  R.  Andrews,  "CTOD  and  J-Relationships  for  a  Growing  Crack",  In  The  Crack 

Tip  Opening  Displacement  in  Elastic-Plastic  Fracture  Mechanics, 

K.  H.  Schwalbe  (ed.),  Springer- Verlag,  Berlin,  p.  191,  (1985). 


surface 


00.'  6  OO'a  00'£  00 ' 9  OO'g  OO'F  OO'E  00'2  00  '  I  OO'O 

(E**OI  X) 

(  N  1  psoi 


uilLL 


Fiqum; 


00 ' 6  00*9  00  'L  00 ' 9  00'S  00 'fr  OO'E  OO'S  00‘T  OO’O 

(E##OT  X) 

[  n  1  peon 


117.2 


203.2 


254  »  10  inch 


Dimensions  of  the  employed  Compact  Tension  specimen  (only  one 
half  of  the  specimen  is  shown,  all  quantities  are  given  in 
millimeter) . 


V  .v  ,V  •- 


Finite  element  modeling  of  the  Compact  Tension  specimen 


Pintle  element  modeling  of  the  Compact  Tension  specimen 


Finite  element  modeling  ot  the  Compact  Tension  specimen 


I 


\  =  0.57735 


\  =.  -0.57735 


T)  =  0.57735 


11  -  0.57735 


Locations  of  the  Gauss  quadrature  integration  points  within  an  eight 
d.  o.  f.  bilinear  isoparametric  element. 


crack  growth 


crack  blunting 


▲a 


Typical  J  •  resistance  cun/e. 


Crack  Growth 


^  original  crack  tip 

v 

current  crack  tip 


24:  "Tangent"  -  definition  of  the  crack  tip  opening  displacement. 


'  tin  *'*- •*-*- J 


(coarse 


crack  growth). 


OHRQlJS  VERSION  <1-5-175 


HUMOUS  VERSION  ^  5-175 


SPECIMEN 

INLUtMtNr  1  HtlHUlJS 


SPECIMEN 


SPECIMEN  FINE  MESH 


Crack  lip  region  (line  mesh) 


7MEN  FINE  MESH 


T-*#OI  X) 

oeidsiQ  6uiuado  d(j_  ^oejQ 


if iT«M»l 


opening  displacement  vs.  crack  growth  (line  mesh). 


Strain  in  y  direction  along  the  line  y  =  0°  for  the  advancing  crack  Radius  r  [  mm  ] 

iarse  mesh). 


Stress  m  y  direction  along  the  line  B  =  0°  for  the  advancing  crack  Radius 


Stress  in  y  •  direction  along  the  line  ©  =  0°  lor  the  advancing  crack 


representation  of  the  von  Mises  equivalent  stress  a 


Log/log  representation  of  the  von  Mises  equivalent  stress  along  line 
w  -  0°  lor  the  advancing  crack  (coarse  mesh). 


advancing  ci 


[siu/m  ]  1  M  ^llsuaQ  A6jaug  uiej;$  DijsBig 


,  «■,  «\ 


.  y  .'■r-'/o";-'.’/ 


Elastic  strain  energy 


Plaslic  strain  energy  density  along  line 


o  _ 
J  o 

,  T3  O 

■  1  ffl  tr 

GC  ’ 


00  ’  6  00 ' 8  00  'L  00 ' 9  00'S  00't>  OO'E  00‘S  00‘t  00‘0 

(I**Ot  X) 


[EUJ/riA|  ]  |dM  M!SU3Q  <6j3U3  UIEJ1S  3llSB|d 


f  iuui t*  ro_  Plastic  stiam  energy  density  along  line 


crack  lip 


HC,  I  MEN 


IMEN 


TablS-L  Relation  between  the  external  (applied)  load  and  the  crack  growth 
(coarse  mesh). 

coarse  mesh  fine  mesh 


F  [N] 

Aa  (mmj 

F(N] 

Aa  [r 

5300 

0 

5300 

0 

6400 

.25 

5520 

.05 

7100 

.5 

5720 

.1 

7650 

.75 

5960 

.15 

8100 

1 . 

6180 

.2 

8250 

1.25 

6400 

.25 

8400 

1 .5 

6540 

.3 

8500 

1.75 

6680 

.35 

8600 

2. 

6820 

.4 

8660 

2.25 

6960 

.45 

8690 

2.5 

7100 

.5 

8710 

2.75 

.  7210 

.55 

8720 

3. 

7320 

.6 

8725 

3.25 

7430 

.65 

8729 

3.5 

7540 

.7 

8732 

3.75 

7650 

.75 

8734 

4. 

7740 

.8 

8732 

4.25 

7830 

.85 

8725 

4.5 

7920 

.9 

8715 

4.75 

8010 

.95 

8700 

5. 

8100 

1 . 

8680 

5.25 

8650 

5.5 

8610 

5.75 

8570 

6. 

8530 

6.25 

8490 

6.5 

8450 

6.75 

8410 

7. 

8370 

7.25 

8330 

7.5 

8290 

7.75 

8250 

8. 

Table  5:  J  -  integral  over  crack  growth. 


J  -  integral 

88.437 

141.117 

193.161 

287.161 
444.633 
536.557 
644.079 
742.951 

847.117 
1008.026 
1133.745 
1251.212 
1366.434 
1469.011 
1554.629 
1614.482 
1620.146 
1624.624 
1642.07 
1618.571 


Iafala  6;  Comparison  of  the  results  for  the  J  -  integral  performed  in  this  work  with 
the  work  of  Hoff  [37]. 


J  at  the  onset  of  crack 
extension. 


slope  of  the  J  -  &a  cun/e 
for  the  first  2,5  mm  crack 
extension. 


this  work 


89.878  N/mm 


393.1  MPa  .using  least 
square  curve  fitting 


200  N/mm 


360  MPa 
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